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INTRODUCTION 


This report is the first annual summary of the research program of the 
Geodynamics Branch. The branch is located within the Earth Survey Appli- 
cations Division of the Applications Directorate of the Goddard Space 
Flight Center. The research activities of the branch staff cover a broad 
spectrum of geoscience disciplines including space geodesy, geopotentiat 
field modeling, tectonophysics, and dynamic oceanography. The NASA pro- 
grams which are supported by the work described in this document include 
the Geodynamics and Ocean Programs, the Crustal Dynamics Project, the 
proposed Ocean Topography Experiment (TOPEX), and Geopotential Research 
Mission (GRM). The reports highlight the investigations conducted by the 
Geodynamics Branch Staff during calendar year 1982. The individual papers 
are grouped into chapters on Crustal Movements, Global Earth Dynamics, 
Gravity Field Model Development, Sea Surface Topography, and Advanced 
Studies. Further information on the activities of the branch or the par- 
ticular research effort described here can be obtained through the Branch 
Office or from the individual staff members. 
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CHAPTER 1 
CRUSTAL MOVEMENTS 


OVERVIEW 


Studies of crustal movements are devoted to the interrelated goals of 
determining tectonic plate motions,, understanding their driving forces 
and the nature of plate interactions, and predicting the occurrence of 
earthquakes. The analysis of global and regional scale crustal 
movements provide information on seafloor spreading and continental 
drift, the creation and destruction of crustal material, and the 
structure, rheology, and dynamics of the earth's interior. Crustal 
deformation in seismically active zones provide information on strain 
accumulation and release, the state of stress in the earth's crust, and 
the mechanisms responsible for earthquakes. The work of the Geodynamics 
Branch includes activities in both the Crustal Dynamics Project which 
involves space geodetic observations and in the Geodynamics Research 
Program which involves basic studies into the dynamic processes of the 
earth. 

This section of the Geodynamics Branch report includes papers on the 
measurement of crustal movements, modeling of earthquake related crustal 
deformations, modeling of sublithospheric stress fields, and the 
development of satellite tracking facilities for space geodetic 
measurements. The four reports are summarized as follows: 1) The 
first report is devoted to observations of crustal movements using 
satellite laser ranging measurements to the Lageos satellite. Repeated 
determination of the distances between observing sites on the North 
American, Australian, Pacific and South American Plates are being used 
to determine both relative plate motion and intraplate rigidity. 2 ) The 
second report focuses on the modeling of crustal deformation subsequent 
to an earthquake. The analyses is based on the physics of stress 
redistribution with time and distance from the seismic event and yields 
information on the rheology and structure of the earth's interior. 3) 
In the third report, the sublithospheric stress field in the Raikal- 
Stanovoy region of Siberia is derived from an analysis of the satellite 
determined gravity field. The computed stress field is interpreted in 
terms of the regional seismicity, and geologic structure and the 
phenomenon of gravity sliding of seismotectonic blocks. 4) The fourth 
report describes data confirming the successful conversion of the lunar 
laser ranging system of the University of Hawaii at Maui into a system 
capable of tracking near-earth satellites. The site is being used in 
measurements of tectonic plate motion involving the Pacific Plate. 

Contributors to this section are: n.C, Christodoulidis, S.C. Cohen, 
R. Kolenkiewicz, H.S. Liu and D.E. Smith. 
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GLOBAL GEOTECTONICS FROM SATELLITE LASER RANGING 

Demosthenes C. Christodoulidis 
and 

David E. Smith 


OBJECTIVE 

The detection of global geodynamics and tectonic processes through the use 
of satellite laser ranging data. 

BACKGROUND 

The accuracy with which tectonic plate velocity vectors can be computed 
increases with the , repetition rate of the measurements and accuracy in- the 
measurement of interplate distances and heights. As a result, a fair amount 
of GSFC activity has been devoted to increasing the temporal resolution and 
precision of baseline and height determinations for laser instruments 
deployed surrounding tectonic plate boundaries. 

RECENT ACCOMPLISHMENTS 

Dynamic reduction , of Satellite Laser Ranging (SLR) provides a unique tool 
for the precise positioning of a point on the Earth's surface in three dimen- 
sions* The sensitivity of the SLR data requires that, for a well conditioned 
solution, the longitude of one of the tracking stations be constrained to an 
a priori value. The vertical components of the stations, however, through 
the tracking geometry provided by the global network and the accurate know- 
ledge of orbital dynamics, are uniquely related to the center-of-mass of the 
earth.*.' .... . . . 

The Lageos laser data for the period May 1976 to April 1982 have been ana- 
lyzed for the determination of sets of annual station locations and earth 
orientation and force model parameters. The results of this analysis, refer- 
red to as SL5, are approaching the useful accuracies required by geophysics. 
The annual height. variations for eight base stations, spread over four tec- 
tonic plates are shown in Fig. 1. Station positions determined from data 
collected prior to 1979 are of inferior accuracy due to the poorer quality, 
quantity and distribution of the data . The three sigma formal errors of 
these heights are between 2 and. 5 cm. Inspection of the heights versus 
their sigmas shows that, for most of the cases, the variations are insigni- 
ficant, with the exceptions being Maui, Hawaii and Yarragadee, Australia. 
The horizontal positions of the stations have an accuracy of better than 2 
cm which has permitted preliminary comparisons of SLR tectonic plate motions 
with those estimated by .Minster and Jordan (1974, 19/8). These comparisons 
are discussed below. 
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SIGNIFICANCE 

The change in the SLR chord lengths, free -rom height variations give 
the horizontal component of the velocity vectors between the tectonic 
plates. These ellipsoidal chords interconnect the projections of the 
tracking station positions onto the reference ellipsoid. As such, they 
may be compared directly with the horizontal changes suggested by Minster 
and Jordan. The Lageos data analyzed was collected from 1979 through 1981. 
In order to increase the sample size, the chords were averged per tectonic 
plate. The computed inter-plate' deformations are shown in Figure 2. The 
heavy lines connect the geographic barycenter of the sites where more than 
one site was available per plate. For the purpose of comparisons the same 
per plate averaging was applied to the Minster and Jordan rates which are 
Shown in Figure 2 also. The agreement between the two independent sets is 
within one sigma of the laser data (approximately + 2 cm/yr) the exception 
being the Australia to S. American rate. It has to be kept in mind, however, 
that this chord crosses two plates (the Nazca plate and the Pacific plate) 
and therefore what we see is the result of the combination of movements 
from four plates. Furthermore, the Minster and Jordan rates reflect 
average movements over mi 11 ions -of -years whereas the SLR rates are 2 year 
averages. Although these results seem to be very encouraging indeed, much 
attention will be given to future higher quality laser results. These 
results will increase the length of the movement records for the laser 
sites and will thus enable the verification of the observed deformation 
rates. 

FUTURE EMPHASIS ■ • 

The results discused above are viewed as preliminary. Studies are now 
underway to better understand and evaluate the variation of the .heights 
and chord lengths and to further improve their temporal resolution and 
accuracy over the entire history of the Lageos mission. Improvements in 
the current data reduction approach include the incorporation of the 
Wahr nutation series, the development and use of a Lageos tidal model and 
Lageos* tidal displacement Love numbers, and the description of the site 
displacement due to ocean loading via an available model. Studies are 
also underway to better model and understand the along-track acceleration 
of the spacecraft, as well as forces due to the earth's albedo. Preli- 
minary error analyses indicate that the incorporation of improvements 
along with better control of the overall systems performance could decrease 
the average positional uncertainty by about 25 percent. 

REFERENCES 

Minster, J.B., T.H. Jordan, P. Molnar, and E. Haines, "Numerical Modeling 
of Instantaneous Plate Tectonics," Geophys. J. of the Royal Astron. Soc., 
36, 541-576, 1974. ; 1 v 

Minster, J.B. and T.H. Jordan, "Present-Day Plate Motions, J. Geophys. 
Res., 83 (Bll), 5331-5354, 1978. 
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CRUSTAL DEFORMATION DURING THE EARTHQUAKE CYCLE 


Steven C. Cohen 


OBJECTIVES 

• . i ! 

The purpose of this work is to develop an understanding of how the earth's 
surface deforms during the cycle of stress accumulation and release along 
major faults * The major research thrusts are the development of numerical 
models of the spatial and temporal patterns of crustal displacements and 
strains and the determination of numerical values of the geophysical para- 
meters that describe the earth's rheology and structure. 

BACKGROUND 

Vertical And horizontal crustal deformations are common characteristics of 
the earth in zones of active seismicity. Certain temporal and spatial pat- 
terns of crustal deformation evolve as the earth passes through successive 
stages of the earthquake cycle: strain accumulation, preseismic movement, 

earthquake rupture, and postseismic rebound. The patterns that develop 
depend on the interaction of the various forces responsible for driving 
earthquakes and tectonic plate movements, the rheology and structure of 
the earth, and the fault geometry. The research reported herein involves 
the development of numerical models that replicate the patterns of defor- 
mation of surface near- and far-field features. The computer Simula* 
tions of an earthquake cycle are fit to observed geodetic data to deter- 
mine key physical parameters of the earth. The specific models that are 
discussed in this report involve postseismic deformation due to visco- 
elastic flow in subcrustal layers of the earth. The flow is activated by 
the coseismic stress field established by an earthquake. As the visco- 
elastic creep advances in the subcrustal layers coupling to the crust 
results in surface deformation. 

RECENT ACCOMPLISHMENTS AND SIGNIFICANCE 

A finite element computer code has been used to study postseismic rebound on 
strike-slip and dip-slip faults. A multilayered viscoelastic rheological 
structure for the ea, th has been used. In this model the earth's lithosphere 
is divided into two layers. The upper layer is elastic and contains within 
it a fault that pentrates through all or only a part of the layer. Below 
this layer is a viscoelastic lower lithosphere represented by a standard 
linear solid rheology. This representation of the lower lithosphere per- 
mits modeling the following deviatoric behavior: elastic deformation 

coseismic to the earthquake fault slip, viscoelastic rebound subsequent to 
the earthquake, and permanent elastic support of stress with a shear 
modulus reduced from the instantaneous coseismic value. Below the lithos- 
phere lies an asthenospheric layer represented by a Maxwell viscoelastic 


1-6 


original page S 

OF POO£ QUALITY 


fluid. The relaxation time of the Maxwell fluid is significantly less 
than that of the overlying linear solid due to the expected thermal 
enhancement of creep processes at depth. 

For strike-slip faults the model predicts that postseismic deformations 
increase as the ratio of fault depth to elastic lithosphere thickness 
increases. The effects of partial relaxation of the lower lithosphere 
are to increase the magnitude of the postseismic displacements and 
strains and to extend the duration of postseismic straining. In addi- 
tion, the graded decreasing viscosity profile with depth retards the 
lateral growth of the zone of significant postseismic restressing. In 
general the effects of asthenosphere flow are significant if the 
asthenosphere depth is less than about 3 fault depths from its surface 
and the viscosity is on the order of 10^0 poise. The partial 
relaxation of the lower lithosphere is significant if the relaxation 
time is not much more than an order of magnitude greater than that of 
the asthenosphere and if the fractional change in the effective 
shear modulus with time is greater than about 1/4. 

The dominant effect in the postseismic rebound from dip-slip event is 
subsidence centered over the fault and uplift in the intermediate or 
far-fi eld region. Over a wide range of dip angles . the magnitude of 
the. subsidence is greatest when the fault extends through one-half the 
Ifthctphere. For greater or lesser fault depths the subsidence is 
r^ifpLed; it becomes uplift when nearly all the lithosphere is pene- 

trated. The width of the subsidence zone is controlled primarily by 
the depth of the top of the asthenosphere and secondarily by the depth 
of the lower lithosphere. The influence of relaxation of the lower 

lithosphere on the surface deformation is much less than that for the 
strike-slip case. Rather the deformation due to asthenosphere relaxa- 
tion dominates the early postseismic vertical displacement pattern 
and continues to develop over times long compared to the Maxwell time. 
This viscoelastic model has been compared to geodetic data obtain in 
1900, 1934, and 1974 following the 1896 Riku-u, Japan earthquake (see 
Figure 1). The deduced value of total lithosphere thickness of about 
40 km and asthenosphere viscosity of 5 x lO 1 ^ poise are moderately 
constrained. The thickness and viscosity of the lower lithosphere of 

20 km and 1 x 1021 poise respectively fit the data but are poorly 

constrai ned. 
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FUTURE EMPHASIS 

The modeling effort will be extended in the future to include strain 
accumulation during the interseismic phase of the earthquake cycle. 
It will be possible to model both the deformation due to a single 
earthquake cycle or to a sequence of recurring events. The finite 
element computer code will provide a general purpose modeling capa- 
bility that will permit undertaking numerical experiments using dif- 
ferent physical assumptions. The deformations that are due to lateral 
and vertical variations in rheological structure will be deduced. 
Numerical experiments will be performed to test the resolvability of 
different models of crustal deformation. The analysis will focus on 
how the surface deformations are influenced by: 1. a thin channel 
asthenosphere; 2. a subsurface low viscosity inclusion; 3. lithos- 
pheric thinning in a fault zone, and 4. elastic creep. 

REFERENCES 

Cohen, S.C., "A Multilayer Model of Time Dependent Deformation Fol- 
lowing an Earthquake on a Strike-Slip Fault," J. Geophys. Res., 87, 
5404-5421, 1982. — ~ 

Cohen, S.C., "Postseismic Deformation Due to Subcrustal Viscoelastic 
Relaxation Following Dip-Slip Earthquakes," GSFC TM 84005, September 
1982; submitted for publication. 

Cohen, S.C., "Discussion of Models of Postseismic Rebound," ESS , 63, 
429, 1982. ~ 
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GRAVITY SLIDING: CRUSTAL DEFORMATION AND SEISMOTECTONIC 

BLOCK MOVEMENTS IN SIBERIA 

Han-Shou Liu 


OBJECTIVE 

Geodynanrics seeks to ascertain the causes of crustal deformation and its 
fundamental connection with geophysical and seismological features. The 
objective of this study is to develop a genetic model which may account 
for the Baikal rift system and the associated seismotectonic block 
movements of Stanovoy Range in Siberia. 

BACKGROUND 

The knowledge of the seismic activity and tectonic structure of the 
Baikal -Stanovoy region in Siberia obtained by the Russian geoscientists 
is profound. It provides geological and geophysical facts to test 
whether the solutions of the focal mechanisms for earthquakes, the 
distribution of orogenic activity,, recent crustal block movements and 
upper mantle structure in this region can be integrated within a single 
Space Geodynamics Program. 

The program for satellite measurement; of the earth's gravity offers a 
variety of models to interpret the interior of the earth. Among the 
various interpretations, the sub-lithospheric stress field serves as 
an intraplate character which corresponds to the state of gravity 
sliding. In this study, the analysis of the Goddard Gravity Field 
Model is used in determining the horizontal projections of lithospheric 
sliding under the Baikal -Stanovoy region. This study involves a search 
for mutual agreement over large areas of the geopotential field with 
respect to seismic stresses. 

RECENT ACCOMPLISHMENTS 

The gravity sliding forces can be inferred from satellite gravity data. 

We have obtained stresses at 1,000 grid points in Siberia caused by 
gravity sliding. The results have determined a seismic gap and recog- 
nized two types of stress conditions in the structure of the Baikal - 
Stanovoy seismic belt. The method developed for identifying stress 
zones from satellite gravity data makes it possible to direct the 
search for the mechanisms of earthquakes and rifts toward these stress 
structures. An interdisciplinary investigation of the stress pattern 
reveals an anomalous mantle beneath the Baikal Rift System. 
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SIGNIFICANCE 


The stresses as shown in Figure 1 are comprised of several consistent stress 
patterns including the remarkable t^nsional stress field under the Baikal rift 
zone and the North-South compression in the Stanovoy Range region. By analyzing 
the satellite stress fields in conjunction with tectonic and seismological 
structures in Siberia, we have obtained the following significant results: 

1. Seismicity and seismic gap: All the strong earthquakes in Siberia are 

concentrated within the Baikal -Stanovoy Seismic Belt (Figure 2) which includes 
the Baikal rift zone and extends eastward as far as the sea of Okhotsk. Around 
118°-120°E there is a seismic gap where earthquake foci are absent. The western 
segment coincides with the Baikal rift zone in which strucures are under exten- 
sion. The eastern segment of the belt is located between 120°E and the Okhotsk 
Sea coast. Within its limites, Baikal -type grabens disappear completely and 
the uplifted topography here points to compression rather than extension. The 
stress pattern, as inferred from satellite gravity data, is in good agreement 
with these seismological features. The most likely point of no strain along the 
Baikal -Stanovoy seismic belt is calculated to lie at 56°N 116°E (Figure I), 
which lies at the position of the aforementioned seismic "gap." 

2. The computed point of no strain (56°N 116°E), or seismic gap, coincides 
almost exactly with the position of the pole of rotation between the Eurasian 
and Amurian plates. This result suggests that gravity sliding may be the main 
cause for seismotectonic block movements. 

3. Figure 1 shows an anomalous lens of upwelling mantle rocks under the rift 
zone which may be verified by seismic or magnetotel luric probing. 

4. Mechanism for Baikal rift system - Gravity sliding in this region, as infer- 

red from satellie gravity data, has provided a mechanism for continental rifting. 
It has been shown that rift evolution in continents is driven by the inexorable 
dissipation of thermal energy from the mantle and the concomitant upward tran- 
sport of the asthenosphere. It may be probable to state that the crustal 

tension causing the Baikal rifting was originated by either the suction force 
acting on overriding plates at trenches or the collision force between the Indian 
and Eurasian plate along the Himalayan frontal thrust. 

FUTURE EMPHASIS 

Satellite capabilites for gravity measurements have been extended so rapidly in 
the past few years that there has not been time to bring them fully to bear on 
the interesting fundamental problems for which they seem so well suited. By 

applying the same method, the author plans to investigate gravity sliding in 
Europe. The gravity sliding stresses in the European crust will be obtained 

form satellite gravity data. The stress calcualtions will provide a real chal- 
lenge to the best geophysical theories, and seismological observations. 

REFERENCES 

Liu, H.S., 1982. Geodynamics of the Baikal -Stanovoy Seismic Belt, Physics of the 
Earth and Planetary Interiors, Vol . 31, 77-82. 
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Figure 1. Subcrustal stresses beneath the Baikal-Stanovoy region as inferred from satellite gravity data. The stress 
pattern is consistent with the anomalous lens of upwelling mantle rocks under the Baikal rift region and the 
associated crustal compression of the S t anovoy Range. 



90 100 110 120 130 140 

Figure 2. Seismicity and stress orientation! of earthquake foci along the Baikal-Stanovoy seismic belt. 
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THE VERIFICATION OF THE LURE SYSTEM AS A SATELLITE LASER RANGING STATION 


Ronald Kolenkiewicz 


OBJECTIVE 

The objective is to qualify the LURE (Lunar Ranging Experiment) laser ranging 
system of the University of Hawaii located on Maui, Hawaii to verify that it 
is suitable to be used as a near-earth satellite system. This enables it to 
join the NASA laser ranging network and thus contribute to the Crustal 
Dynamics Project. 

BACKGROUND 

In the past the LURE system was designed exclusively for laser ranging to the 

moon. NASA was in need of a Satellite Laser Ranging (SLR) system on the 

Pacific Plate in order to perform Global Studies for the Crustal Dynamics 
Project. These studies include monitoring of tectonic plate motion as well 
as polar motion and earth rotation. The LURE system was modified to enable 

it to track near-earth satellites. In order to verify the ability of the 

LURE system to accomplish this goal, a SLR system, MOBLAS I was collocated 
in the vicinity of the LURE laser. 

These collocated systems tracked the BE-C, Starlette, and LAGEOS satellites 
to demonstrate that the LURE system gave the same ranging information as 
MOBLAS 1. These tests were conducted during September through November 1981. 

RESULTS 

The laser ranging data from the two collocated laser systems (LURE and MOBLAS 
I) was analyzed by selecting simultaneous passes of data to the same satel- 
lite. The residuals to an orbit calculated to these data based, upon the 
known survey position of the sites, were analyzed for relative range and 
timing biases. Figure 1 shows the laser range residuals i "om collocated 
lasers, about an orbit. The average residual for each station is deter- 
mined. These averages are then differenced to give the mean residual 
difference for that pass. Figure 2 shows these mean residual differences 

plotted for the 36 satellite passes analyzed. Initial results indicated 
there was a relative bias of 5 cm between the stations. Allowing for one 
station to move to a new position relative to the other removed this bias. 

A survey confirmed these analyses by showing the original position of 

MOBLAS I was in error by 6.8 cm in latitude and 3.8 cm in longitude and 

the calibration board distance was in error by 1.1 cm. 

When the corrected position was used in the anlaysis the measurments made by 
both systems were in agreement with each other by less than 1 millimeter. 
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Figure 2 presents the results of this experiment. The mean residual difference 
for the 36 passes analyzed are seen to be spread between plus and minus 4 cm 
with an average of zero and a sigma of 2.4 cm. 

CONCLUSIONS 

The LURE system is capable of satellite laser ranging at the centimeter level, 
and SLR systems can detect survey errors at the centimeter level. This 
implies that relative site positions can be measured at the centimeter level 
providing tectonic and earth orientation measurements of value to the Crustal 
Dynamics Project investigators. 
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CHAPTER 2 


GLOBAL EARTH DYNAMICS 


OVERVIEW 

Studies of global earth dynamics are devoted to improving the understanding 
of the origin of motions of the whole Earth and determining the structure and 
mechanical properties of the Earth's interior. The techniques employed in 
such studies are as varied as the studies themselves. They range from the 
determination or changes in satellite orbits, to spectrum analyses of polar 
motion data, to the inversion of gravity field representations. The reports 
in this chapter summarize as follows five research activities underway in the 
Geodynamics Branch during 1982: 

1) The first report focuses on the determination of the effective viscosity 
of the lower mantle using observed accelerations in the node of Lageos' 
orbit. The results are discussed in terms constraints on the allowed 
modes for mantle convection. 

2) The second report is the first of two papers on polar motion. Data from 
both the International Latitude Service and observations of Lageos have 
been analyzed to determine both the annual and Chandler Wobble compo- 
nent of polar motion. For both of these terms the frequencies, ampli- 
tudes, phases, and Q's of the signal have been determined. 

3) Continuing the study of polar motion, the next report discusses the 
development of a matched filter to improve the signal recovery from a 
noisy time seris. This filter is being used to determine the frequency 
and Q of the Chandler Wobble. 

4) The fourth paper describes the use of gravity data in an information 
theory based inference scheme to determine the lateral density distri- 
bution of the Earth. 

5) Finally, a comparison is presented on length-of-day determinations from 
atmospheric wind data and from Lageos observations. The correlation 
between the 2-day sets argues that changes in the zonal winds are a 
major contributor to length-of-day changes,, 

Contributors to this chapter are B.F. Chao, D.C. Christodoulidis, R. Gross, 

D. Rubincam, D.E. Smith, and George Wyatt. 


2-1 


ORIGINAL FAGS IS 
OF POOR QUALITY 


POSTGLACIAL REBOUND ORSERVED BY LAGEOS AND THE EFFECTIVE VISCOSITY 

OF THE LOWER MANTLE 

David P. Ruhincam 


OBJECTIVE 

The objective of this work is to infer the effective viscosity of the lower 
mantle from observations of the node of Lageos' orbit. 

BACKGROUND 

The Lageos satellite was launched in May 1976 for the purpose of studying 
crustal motion, plate motion, polar motion, and earth rotation (Smith and 
Dunn, 1980). Five years worth of laser ranging data reveal an acceleration 
in the node of Lageos' orbit which is not currently modeled in the GEODYN 
orbit determination computer program (see Fig. 1). Its value is p = -8.1 
± 1.8) x 10“° arcseconds day" 2 ; hence the rate at which the node is 
progressing along the equator (p = 0.343 degrees day -1 ) is becoming smaller. 

An obvious geophysical interpretation of the data is to assume that the 

. _ . * I K > • >(•<• • ... . • 
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RECENT ACCOMPLISHMENTS 

A complete description o^ the work summarized here may be found in Rubincam 
(1983). Briefly, the observed J^ cannot be due to the tidal deceleration of 
the earth and the subsequent relaxation of the earth; it is a factor of 40 
too small. Also, the long-period 18.6 year ocean tide does not appear to be 
the cause, since it woud have to be far from equilibrium. 

The acceleration of the node does appear to be due to the postglacial rebound 
of the earth. During the Quaternary T.ce Age the weight of the ice caps in 
the polar regions squashed the earth into a more elliptical shape; it is now 
relaxing and thus decreasing J 2 . The magnitude of J? is controlled in part 
by the effective viscosity of the mantle. Hence the observed Jg places 
constraints on earth structure. 

Wu and Peltier's (1982) LI and L2 earth models plus an assumed ice cap 
history were used to compute J 2 values to compare with the observed value. 
The LI model has a 120 km-thick lithosphere, a lf)21 Pa s mantle, and an 
inviscid core. The L2 model is the same as the LI, except that the effec- 
tive viscosity of the lower mantle (below 671 km depth) is lf)22 p a s . 
The resulting J 2 values were -9.8 x 10“ 19 s"* 1 for the LI model and -32.0 x 
10-19 S -1 f or the L 2 model. Clearly the LI model is closer to the observed 
value of (-8.1 ± 1.8) x 1D-I q s-1 and is preferred to the L2 model. 


earcn's equatorial ouige is 
tional to the J 2 term in 
gravitational potential, the 
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getting smaller. Since the rate p is propo 
spherical harmonic expansion of the earth 
observed p = (8jl±1.8) x 10- fi arcseconds 0 day 


implies that J 2 is changing at the rate of J ? = (-8.2 ± 1.8) x ID 
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SIGNIFICANCE 


The results indicate that the effective viscosity of the mantle must be a 
fairly uniform 10 2 * Pa s. In other words, there is no great difference 
between the effective viscosities of the upper and lower mantle. Hence any 
objection to who! i mantle convection based on such a difference apparently 
cannot be invoked. 

FUTURE EMPHASIS 

The behavior of the node of Lageos' orbit will be monitored in the future to 
see whether it still conforms to that predicted by postglacial rebound. 
Also* other possible causes for the acceleration of the node will be 
investigated. 
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DYNAMICS OF THE EARTH'S POLAR MOTION 


B. Fong Chao 


OBJECTIVE 

To analyze polar motion data obtained by geodetic techniques including classical 
astrometry and satellite laser ranging, and to compare the observations with the 
atmospheric and seismic excitation functions. 

BACKGROUND 

The Earth's rotational axis does not remain fixed relative to the body of the 
Earth. Instead, the interception of the axis with the surface of the Earth 
(f.e., the pole) traces out a quasi -periodic path about some mean position on a 
scale <10 meters. This motion is known as the polar motion of the Earth. The 
polar motion in principle consists of a number of components arising from various 
dynamical processes. The two major components are the annual wobble with a 
yearly period and the 14-month Chandler Wobble. 

The annual wobble is a forced motion and the variations in atmospheric mass 
distribution have long been identified as the major cause for the annual wobble. 
Other mechanisms, such as ocean circulations, core-mantle coupling, solar winds 
and radiations, are believed to have minor effects on Earth's polar motion. 
Numerically, lacking any information that indicates otherwise, the annual wobble 
has been modeled as a pure sinusoid with a period of 365.25 days. 

The Chandler wobble, on the other hand, s the Earth's free motion corresponding 
to the Eulerian nutation of a rigid body. Its excitation sources have remained 
unidentified despite an ongoing effort over many decades. Two plausible candi- 
dates have been proposed: atmospheric and seismic excitations. However, as to 
date, both seeem to suffer from insufficient energy. Despite this difficulty 
(rightfully or not), the period and the 0 of the Chandler Wobble have been esti- 
mated by means of various numerical techniques; and there has been no agreement 
on the structure of Chandler period(s) (whether it is single, multiple, or even 
time-variable), nor on the value of the Chandler 0 (estimates range from from 
~3Q to ~1000, depending on data and techniques used). 

RECENT ACCOMPLISHMENTS AND SIGNIFICANCE 

Naturally, for our purposes (namely, the numerical modeling of polar motion 
using atmospheric and seismic excitation functions), the essential first 
step i$ to separate the annual and the Chandler wobbles into optimum numerical 
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representations— in the present case, complex harmonic functions of time. 
Thus, we completed a harmonic analysis of the following two sets of polar 
motion data: (1) monthly solutions by ILS (the International Latitude 
Service) using classical astrometry, spanning 80 years (1900-1979), and (2) 
5-day solutions of LAGEOS (the Laser Geodynamics Satellite) using satellite 
laser ranging, spanning 5.7 years (1976.4-1982.1). Tfie method used is the 
"Autoregressive (AR) method" (Chao and Gilbert 1980), devised to estimate 
the complex frequency (frequency and 0) and the complex amplitude (amplitude 
and phase) of each complex harmonic function in a time series. Tables 1 and 
2 give the results for ILS data (Chao, 1983a). Major conclusions include 
(i) the annual wobble is stationary over the long period of time spanned by 
ILS observation, and (it) the ILS data are consistent with a Chandler wobble 
model that consists of four components. Tables 3 and 4 give the LAGEOS 
results (Chao, 1983b). The multiple-component behavior of the Chandler 
wobble, if it exists, cannot be resolved due to the short record length. 
Note that the estimated annual periods are significantly shorter than the 
theoretical value 365.25 days (by 5-10 days) during 1976.4-1982.1. It should 
be pointed out that the observations in Table 3, including the anomalous 
periods, the decay in amplitudes, and the differences in amplitudes and 
phases, simply reflect the yearly behavior of the driving mechanism(s) of 
the annual wobble, and can and should be accounted for by the latter. 

On the theoretical front, the problem of excitation of the Earth's (free) 
Chandler wobble by a transient, indigenous source (as exemplified by a seismic 
event whose duration is generally much shorter than the Chandler period) has 
been revisited (Chao 1983c), with particular emphasis on the compatibility in 
coordinate systems used for observations and for computations. It is found 
that the difference between the two types of coordinate systems is far from 
trivial. However, for most practical purposes, this discrepancy is small 
compared with the amplitude of the polar motion itself. 

FUTURE EMPHASIS 

The long ILS data set is ideal for statistical studies for the excitation of 
polar motion. One such study is a statistical analysis of the relationship 
between the polar motion and the occurrence of major earthquakes. Thus, 
correlation between the longitudes of over 1500 earthquakes (over the same 
period spanned by the ILS observation) and the phases of the pole position at 
the time of the events will be calculated in order to detect the existence of 
any possible cause-and-effect relation. 

The high-density, accurate LAGEOS data, whose annual and Chandler wobbles 
having been separated and analyzed, can now be compared with the atmospheric 
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angular momentum data (also available since 1976, see Rosen & Sal stein 1981). 
In particular, we can compute the polar motion that is caused by the atmos- 
pheric variations and subsequently conduct a wobble-by -wobble comparison. 

This presumably will account for the behavior of the annual wobble presented 
in Table 3, and determine how much of the Chandler wobble can be accounted 
for by the atmosphere and hence give clues to the importance of other exci- 
tation agents. 



Period (days) 

0 

Amplitude (") 

Phase (°) 

ILS-X 

365.20 

-1480 

0.0884 

108 

ILS-Y 

365.10 

1370 

0.0840 

10 


Table 1. AR 

estimates 

for the ILS annual Wobble 



Comp. No 

Period (days) 

n 

Amplitude (”) 

Phase (°) 

I 

406.45 

1840 

0,0267 

216 

II 

426.00 

711 

0.1199 

273 

III 

437.46 

-189 

0.0573 

342 

IV 

452.73 

180 

0.0686 

150 



(a) 



Comp. No. 

Period (days) 

0 

Amplitude (") 

Phase (°) 

I 

406.85 

-386 

0.0170 

132 

II 

426.15 

703 

0.1195 

185 

III 

437.43 

-184 

0.0572 

253 

IV 

452.39 

320 

0. 0569 

54 


(b) 

Table 2. AR estimates for the Chandler wobble from 
(a) IIS-X (b) ILS-Y. 
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Period (days) 

0 

Amplitude (") 

Phase (°) 

LAG-X 

360.0 

13.2 

0.153 

196 

LAG-Y 

354.9 

8.7 

0.182 

80 


Table 3. AR 

estimates for the LAGEOS annual wobble. 



Period (days) 

0 

Amplitude (") 

Phase (°) 

LAG-X 

434.2 

22.5 

0.188 

319 

LAG-Y 

436.4 

12.8 

0.245 

238 


Table 4. AR estimates for the LAGEOS Chandler wobble 
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POLAR MOTION 


Richard S. Gross 


OBJECTIVES 

The objective of this research effort is to study some of the outstanding 
problems in the field of polar motion such as the nature of the excitation 
mechanism of the Chandler wobble, the value of the period and 0 of the 
Chandler wobble, the Markowitz wobble and the long term stability of the 
latitude observation stations. 

BACKGROUND 

In 1765 Euler predicted that any rigid body, including the earth, which is 
rotating about some axis which is not its principal moment of inertia axis 
will experience a wobble as it rotates. In 1891 Chandler was able to observe 
this wobble of the earth and since 1899 it has been observed and studied by 
the International Latitude Service ( ILS) and, more recently, by many other 
organizations. This motion of the rotation pole about the principal moment 
of inertia axis is primarily composed of two components. The first has an 
annual period and is a forced motion of the pole commonly believed to be 
caused by meteorological events. The second is Euler's free nutation at a 
period of approximately 434 days and is now known as the Chandler wobble. 
However, the motion of the rotation pole also exhibits a long term drift in 
its mean position and there is some evidence that a component (called the 
Markowitz wobble after the person who proposed its existence) might exist 
having a period of roughly 24 years. 

RECENT ACCOMPLISHMENTS 

As a graduate student I used a damped least squares technique to obtain the 
polar motion along with station corrections from the homogeneous ILS latitude 
observations in order to obtain a polar motion data set at 19-day time 
intervals. This polar motion time series spans the years 1899-1979 and can 
be used to study the long term behavior of the motion of the pole. Since 
the solution technique took into account the individual station contributions 
via the station corrections we are also able to study the long term stability 
of the individual ILS latitude stations. 
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As our first use of this newly derived polar motion time series we will 
ascertain when the Chandler wobble is undergoing excitation by a technique 
based upon matched filters. A matched filter is the optimal linear filter 
that maximizes the signal to noise ratio in the output at some instant. This 
is desirable if we are trying to detect the onset time of some signal embedded 
in noise. The assumption is that our chances of detecting this signal will 
be increased if the signal to noise ratio is a maximufo. For a time series 
consisting of a known signal s(t) of amplitude A and arrival time t 0 embedded 
in noise n(t): 

x(t) = As(t-t 0 ) + n(t) 


it can be shown (see, e.g., Robinson (1967)) that the desired filter in the 
frequency domain is simply S*(w), the complex conjugate of the Fourier 
transform of the signal s(t): 

Y (to) = X(w)S*H 


where X(co) is the Fourier transform of the input signal x(t) and Y(w) 
is the Fourier transform of the output signal y(t). It can be shown (see 
Gross (1982)) that if the signal s(t) we are trying to detect is damped 
harmonic motion then the output signal y(t) will exhibit a peak at the signal 
arrival time t 0 and it will decay exponentially from this peak value. In the 
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peak in the output will be rather broad. Of course, we would like the peak 

to be as delta like as possible. In fact we can form a narrower peak by using 
the filter: 


Y(w) = X(o>) 


$*(w) 


(IINM II* + ||SM I I 2)1/2 


where N(w) is the Fourier transform of our estimate of the noise in the 
input time series. To illustrate the use of the above filter figure 1 shows 
the x-component of a synthetic time series we generated which consists of 
three signals embedded in noise. Figure 2 shows the amplitude of 
the time series resulting from applying the above filter (with N=0.1) to this 
input time series. As can be seen the output is peaked at the signal arrival 
times which are indicated by the dashed lines. 

We are currently applying this noisy matched filter to various polar motion 
data sets. As an aid in interpreting the results we are also using it on 
many synthetic time series. 

FUTURE EMPHASIS 

We shall continue our investigations into the properties of the excitation 
process of the Chandler wobble as well as study other topics in the field of 
polar motion. In particular we will calculate the period and Q of the Chandler 
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wobble by assuming that it is the result of a white noise process exciting a 
resonant system. We will also study the Markowitz wobble as manifested in 
the data set derived by Gross (1982) and we will study the long term stability 
of the ILS latitude observation stations as indicated by their station 
corrections. 
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INFORMATION THEORY LATERAL DENSITY DISTRIBUTION 
David P. Rubincam 


OBJECTIVE 

The objective of this work is to use information theory to infer the lateral 
density distribution inside the earth, using the external gravity field, 
crustal models, and an assumed radial density distribution as data. 

BACKGROUND 

The external gravity field gives information about the lateral density dis- 
tribution inside the earth. However, it does not given enough information 
to determine the actual density distribution pj(Y,e,X). The lateral 
density distribution must be inferred, not inverted from the data. 

The inference technique used here is information theory better known as the 
Maximum Entropy Method (MEM). The rationale for using information theory is 
that it provides, in a certain sense the least prejudiced inference from the 
available data. In other words, it does not read more into the data than is 
actually there. 

RECENT ACCOMPLISHMENTS 

A complete description of the work summarized here may be found in Rubincam 
(1982). Briefly, a density distribution pj(y,9,X) is inferred from the 
the GEM 1DB gravity field coefficients C £m , S £m (complete to degree and 
and order 36) after subtracting off the coefficients of the hydrostatic bulge 
and an assumed crustal model. Two different models were used: one an isosta- 
tically compensated crust, and the other an uncompensated crust. Also the 
radial part of pi(y> 0,*) is taken to be the Parametric Earth Model of 
Dziewonski et al . (1975). 

The following results were obtained. The density anomalies generated by 
pi(y,Q,x) tend to be spread throughout the mantle, but the largest 
anomalies occur near the top of the mantle. Typical density variations in 
the mantle are +0.D04 g cm “ 3 when the compensated crust is removed and 
+0.04 g cm -3 when the uncompensated crust is removed. The lateral density 
variations across the equatorial plane are shown in Figs. 1 and 2 for compen- 
sated and uncompensated crusts remoived, respectively. 

SIGNIFICANCE 

The pj(y,0,X) for either compensated or uncompensated crusts shows no 
obvious convection pattern in the mantle (such as organized convection cells). 
Hence mantle convection cannot be inferred from the gravity field data alone. 
Also, the information theory density distributions are quite different from 
the minimum elastic energy solution of Kaula (1963) and from the minimum 
elastic plus gravitational potential energy solution of Sanchez (198D). 
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FUTURE EMPHASIS 

The extension of the information theory approach so as to include seismic 
data appears quite difficult. Hence no future work in this area is planned. 
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Figure 2. Lateral density variation in the equatorial plane of the earth for the case of 30 km of uncompensated 
crust removed. Sea level is displaced from the eartYs surface for ciarity, and the rock-equivalent 
topography is g>eatly exaggerated. Anomalies in the core are not shown. 
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ATMOSPHERICALLY INDUCED VARIATIONS IN THE LENGTH-OF-DAY DERIVED FROM LAGEOS 


David E. Smith 

Demosthenes C. Christodoulidis 
George H. Wyatt 


OBJECTIVE 

A comparison is made of the length-of-day derived from Lageos with the length- 
of-day inferred from atmospheric wind data in order to determine the extent of 
the correlation, accuracies of the two data sets, and to better understand 
the causes of the changes in the earth's rotation rate. 

BACKGROUND 


Recently, global atmospheric wind data have been used to derive estimates for 
the angular momentum of the atmosphere (Rosen & Salstein, 1981) and, under the 
assumption of the conservation of angular momentum between the solid-earth and 
atmosphere, the resulting change in !’ ■*» rotation rate of the solid earth has 
been derived. In these computations tne atmospheric zonal winds have been 
integrated from the surface (1000 mbars) to the 100 mbar level (i.e. 16 km) and 
averaged into latitude bands from which the angular momentum of the atmosphere 
has been derived every 12 hours since 1976. Considerable seasonal and short 
period (tens-of-days) variations are evident in this data set which has been 
compared with the actual changes in the earth's length-of-day derived from 
laser tracking of the Lageos spacecraft between May 1976 and July 1981. 

RECENT ACCOMPLISHMENTS 

The atmospheric angular momentum data was averaged to 5-day values so as to 
coincide in time with the Lageos length-of-day (L0D) values. The two data sets 
were analyzed using identical methods. A spectral analysis of the atmospheri- 
cally inferred L0D for the whole 6-year period showed considerable power at a 
period of 1 year and a smaller, but significant, amount at 6 months. There 
appeared to be no secular trend or very long period evident in the data. 
The spectrum of the Lageos L0D showed greatest power at 6 months and almost 
the same at 1 year and significant power also at about 1166 days and 285 
days. These latter periods are periods related to the tidal perturbations 
of the orbit by the and S 2 tides respectively and are not actual var- 
iations in the rotation rate of the earth. They were estimated and removed 
from the Lageos data set. Tne Lageos data set also showed a secular trend 
that indicated that the length-of-day was decreasing (earth was speeding up) 
at a rate of about 0.15 msecs/yr over the 6-year period. This trend is 
believed to be a result of changes within the earth's core and probably 
t>v,iporary in nature. Historical records of the earth's rotation rate show 
"decade length" changes that have generally been ascribed to core motions. 
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After removal of the secular trend and the tidal perturbations of the Lageos 
orbit (1166 days, 285 days) from the Lageos LOD data, the two data sets were 
compared. Figure 1 shows the atmospheric data superimposed on the Lageos data 
after the average value has been removed from each set* The correlation be- 
tween the two "raw" data sets is 0.82 and the rms difference is 0.19 milli- 
seconds. A major contributor to the rms difference is the very large noise 
present in the 1976 and 1977 Lageos LOD values. The best year, in terms of 
Lageos tracking, is 1980 when the correlation between the two LOD data sets 
is 0.90 and the rms difference in LOD values is 0.15 milliseconds. 

The high degree of correlation between the two data sets is partly due to the 
strong annual and semi-annual periods in both data sets, and the rms difference 
between the curves is also greatly affected by the fact that each data set has 
different amplitudes for these terms. When the annual and semi-annual terns are 
removed from the data for the year 1980 the correlation falls to 0.67 (stil : 
very good) and the rms difference falls to 0.10 mil liseconds. The decrease in 
rms can be almost completely explained by the removal of the power at 1 year 
and 1/2 year that appears in the single point rrns. The true level of agreement 
between the atmospheric and Lageos data sets can be best seen in the year 1980 
when a comparison of the two data sets is made after they have been weakly 
smoothed (so as to remove some of the measurement noise) and the anpual and 
semi-annual terms removed. This comparison is shown in Fig. 2. The smoothing 
is achieved by a Gaussian weighting where the FWHM is about 8 days. The 
resulting rms difference is 0.069 milliseconds and the correlation is raised to 
0.78. 

CONCLUSIONS AND FUTURE EMPHASIS 

Comparisons of the atmospheric and Lageos data sets strongly support the con- 
cept that most changes in the earth's length-of-day can be ascribed to 
changes in the earth's zonal winds. Good correlation appears to exist be- 
tween the data sets for periods of about 10 days to more than a year. The 
disparity that appears in the amplitudes (0.17 msecs) at 1/2 year period is 
probably caused by winds not included in the atmospheric data, i.e., those 
above 100 mbars (Lambeck, 1980). At 1 year, the disparity (0.12 msecs) in 
amplitudes is probably associated with the lack of atmospheric data above 
100 mbar but- is also affected by the lack of consideration of variations in 
the angular momentum of the oceans (Eubanks, 1983). 

Future work will try to explain the present differences and investigate other 
causes of change in the earth's length-of-day. 
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Figure 1. Comparison of Atmospheric Angular Momentum with Lageos Variations in LOD. 
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CHAPTER 3 


GRAVITY FIELD MODEL DEVELOPMENT 


OVERVIEW 


Knowledge of the earth's cravity field is fundamental to understanding 
the dynamics of the earth. For solid earth physics, knowledge of the 
variations in the gravity field provides information on the earth's 
physical properties and geodynamic processes and places constraints on 
the internal structure of the earth. In oceanography, knowledge of 
departure of the actual sea surface from a unique equipotential surface 
of the earth's gravity field (the geoid) can reveal information on oceanic 
circulation. In addition, other areas which benefit from knowledge of 
the earth's gravity field are satellite orbit determination and classical 
geodesy. 

This chapter describes: (1) the development of a refined gravity field 

model for LAGEOS in order to obtain more accurate orbits that will 
provide better baseline measurements for the analysis of tectonic plate 
motion; (2) the study of the evolution of the orbit of LAGEOS for the 
purpose of improving our knowledge of the force model for the satellite's 
motion; (3) the development of specialized earth gravity field model of 
the STARLETTE satellite to improve the precision of that satellite's 
orbit, and (4) current and future geodynamic computational systems 
development. 

Contributors to this chaper include D.C. Christodoulidis, F.J. Lerch, 
J.G. Marsh, B.H. Putney and D.E. Smith 
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GRAVITY MODEL REFINEMENT FROM LAGEOS 


Francis J. Lerch 


OBJECTIVE 

The objective is to improve the gravity field model for Lageos in order 
to obtain more accurate orbits that will provide better baselines for 
plate tectonics, improved polar motion and earth rotation. Also the 
purpose is to improve the low degree and order geopotential which will 
provide for refinement of our knowledge of the broad features of ocean 
ci rculation. 

BACKGROUND 

The development of the gravity field for Lageos is a principal investiga- 
tion in the Crustal Dynamics Project. Based upon the previous Goddard 
Earth Model, GEM 9, baselines from Lageos orbits gave an rms of approxi- 
mately ±5 cm from a global set of laser tracking sites. Statistical 
error analysis based upon the GEM 9 accuracy of the coefficients confirms 
that the baseline error (±5 cm) from Lageos can be attributed signifi- 
cantly to the gravity field model. From the accuracy of the Lageos 
laser data (down to +10 cm) error analysis shows that the gravity 
model can be improved on Lageos to the extent that baseline errors should 
reduce by at least 50 percent of previous values when considering errors 
only due to the gravity field. 

RECENT ACCOMPLISHMENTS 

A refined gravity field model, Goddard Earth Model GEM-L2 (Lerch, et 
al.), has been derived using the Lageos orbital data yielding better 
baseline measurements for the analysis of tectonic plate motion. This 
field also contributes to a better understanding of the broad features 
of ocean circulation through its improvement of the long wavelength geoid 
(terms through degree and order 4) with an accuracy estimated at ±8 cm. 
In GEM-L2 two and a half years of Lageos laser data acquired from over 
20 well -distributed stations were combined with the existing data from 
the best satellite-derived model, GEM-9 (Lerch, et al.). Testing shows 
that the Lageos gravity field error at long wavelengths is less than half 
that for GEM 9. A comparison of global laser "base" stations from inde- 
pendent data sets of alternating 15-day data segments over the 2 1/2 
years of Lageos show total interstation position to ±1.8 cm (Table 1) 
when using this new model. The same comparison using the 1979 versus the 
1980 Lageos data yields ±5.2 cm; this difference in agreement re- 
flects the change in data distrbituion and other systematic errors along 
with the tectonic motion which has occurred between these chronologi- 
cally distinct data sets. Polar motion values were adjusted in GEM-L2 
to a precision of 10 cm, which is shown to be necessary to achieve the 
accurate baseline and geopotential results. 
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The accuracy of the GEM-L2 long wavelength geoid of ±8 cm should 
provide improved accuracy in the broad features of ocean circulation from 
application of Seasat altimetry. 

Results based upon the gravity model GEM-12 including adjustments of BIH 
polar motion and earth rotation from Lageos laser data have shown that 
laser tracking systems yield an intrinsic accuracy of 1.8 cm for a global 
set of baselines (Table 1). Tests in this table clearly show that it is 
equally important to adjust the polar motion and earth rotation as well 
as the gravity field with the Lageos data in order to achieve the 1.8 cm 
baseline agreement. 

FUTURE EMPHASIS 

Future work is continuing with use of more recent and accurate laser data 
from Lageos. 
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LAGEOS ORBIT EVOLUTION 


Demosthenes C. Christodoulidis 
and 

David E. Smith 


OBJECTIVE 

Analysis of the evolution of the orbit of the Lageos satellite can improve 
our knowledge of the force model for the satellite's motion. More accurate 
representations of the earth, ocean and atmospheric tides and of the 
gravity model are of general geophysical interest. Improvements in the 
perturbation model for drag or radiation pressure will allow better 
determination of the geodynamic measurements which are the main goals of 
the Lageos analysis. 


n a au/m\ai urn 

dAu^jkOUinu 

Rubincam (1980) has suggested charged particle drag as the cause of an 
observed downward trend in the satellite's semimajor axis. Anselmo, et 
al . (1983) have proposed a model for the variations about this trend based 
on the effect of reflected radiation from an earth whose albedo exhibits 
hemispherical assymetry. 

The procedure employed to determine a set of locations for the network 
of stations over the total time period for which laser t servations of 
the Lageos satellite are available has been described by .alth, et al. 
(1978). It includes the adjustment of Lageos elements at the epoch of 
each 30-day time span over which a dynamically consistent orbit is deter- 
mined, and which is the subject of the orbit evolution studies. 


RECENT ACCOMPLISHMENTS 

The estimated values of along-track acceleration at 30-day intervals 
obtained in a recent solution are shown in Figure 1. Their average effect 
of -3.3 x 10" 12 m/sec 2 corresponds to a semi major axis decrement of 
approxiamtely 1.1 mm/day. The semimajor axis time history due to this 
unmodeled force is shown in Figure 2 as two time spans during which the 
semi major axis decrement was about 1.1 mm/day separated by a period in 
1979 when it was somewhat larger in magnitude. 
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The incl incation and node residual signatures shown in Figures 3a and 3b 
include the effects of errors in the adopted tidal model on the estimated 
orbit parameters. The variations can be reduced by correcting the 
inclination and node values with sinusoidal functions at the three tidal 
frequencies of 1051, 521, and 280 days corresponding respectively to the 
Ki, K2 and S2 tides. If we assume that the earth tidal model {Wahr, 1979) 
is correct, these functions can be used to correct the second and fourth 
degree terms of the ocean tidal model (due to Schwiderski ). The values 
derived in the Lageos analysis are shown in Figure 4 together with those 
corresponding to the Schwiderski ocean tidal model and that due to Yoder 
(JPL E.M. 335-30). 

SIGNIFICANCE 

The relatively low amplitude in the second degree S2 tide observed by 
Lageos could be caused by the atmospheric tide to be expected at this 
frequency. The only other serious disagreement with the ocean tidal models 
is in the phase of the fourth degree term at the Kl frequency. This is a 
manifestation of the long period variation observed in the nodal signature 
of Figure 3b. An error in our model of the 18.6 year tide could provide 
an alternative reason for this error. Although the effect of the neglected 
ocean tide at this frequency is too small by a factor of 10, an increase 
in the effective Earth tidal Love number k? from 0.299 to 0.6 could 
produce a nodal variation similar to that observed. 

The amplitude of a quadratic function of time which can be used to describe 
an observed curvature in the nodal signature corresponds to. a variation in 
the second zonal harmonic amounting to -.35 x 10 -/yr. Rubincam 

(1983) has proposed post-glacial rebound as the cause of this variation, 
which when removed leaves us with a signature dominated by errors in the 
adopted time reference system of the BIH Circular D. 

FUTURE EMPHASIS 

Numerical models of the earth's albedo will be calibrated and employed to 
help explain the observed variations in the semi major axis and eccentri- 
city of Lageos. As the time history of the inclination and node signa- 
ture extends, better separation of the tidal frequencies will be possible. 
Improved definition of the quadratic nodal signature which can be described 
by a variation in the earth's second zonal harmonic can also be expected. 
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Figure 3. Evolution of Inclination and Node Errors 
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Figure 4. Comparison of Ocean Tidal Models 
























GEODYNAMICS AND GEODETIC PARAMETER ESTIMATION 
FROM STARLETTE LASER TRACKING DATA 

Janies G. Marsh 


OBJECTIVES 


Laser tracking data covering a 4-year period (1975-1978) have been ana- 
lyzed for the development of a specialized earth gravity model for 
STARLETTE precision orbit computations. Subsequent to this development, 
a detailed analysis of the STARLETTE laser data recorded during 1980 has 
been performed for the estimation of tracking station coordinates, polar 
motion, GM, and ocean tidal constituents. This solution is a simul- 
taneous least squares solution based upon normal equations for 65 5-day 
arcs of laser tracking data. 


BACKGROUND 


The Starlette satellite launched in February 1975 (CNES/GRGS, 1975) by 
the French Centre National d 1 Etudes Spatiales (CNES) was designed to 
minimize the effects of non-gravitational forces and to obtain the highest 
possible accuracy for User range measurements. The satellite is a sphere 
with a radius of 12 cm. It is a passive satellite constructed of a core 
which is largely Uranium 238 covered by an aluminum alloy shell which 
contains 60 laser corner reflectors. The mass of the satellite is 47.295 
kg. The orbit has perigee and apogee heights of 810 km and 1105 km 
respectively and an inclination of 49.8 degrees. 


Laser tracking data recorded during the time span from 1975 to 1981 have 
been used in these analyses. The data analyzed have been contributed by 
the CNES, Smithsonian Astrophysical Observatory (SAO), Goddard Space 
Flight Center (GSFC) and the Technische Hogeschool Delft in the 
Netherlands. 

Initial orbit computations using 5-day arc lengths of the laser data 
produced rain fits for a typical set of orbits in 1979 which ranged from 
1.5 to 3 meters. In these reductions the GEM 10B gravity model was used. 
Analyses of these solutions indicated that the primary error source was 
the earth gravity model, and thus a comprehensive effort was undertaken 
to adjust all gravity model coefficients which would have an effect of 1 
cm or larger on the STARLETTE orbit. The final set of coefficients 
adjusted consisted of a complete (36,36) model plus selected higher order 
coefficients to degree (48). 
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The high level of orbital accuracy provided by this gravity model solution 
has provided a basis for the adjustment of a global set of tracking 
station coordinates, polar motion, GM, and ocean tidal coefficients. 
Laser tracking data recorded during 1980 have been used in these analyses. 

RECENT ACCOMPLISHMENTS 

The new STARLETTE gravity model designated PGS~1331 reduced the rms 
orbital fits for the laser data (5-day arcs) from several meters to a few 
decimeters. 

The analyses of the more precise globally distributed laser tracking data 
recorded during 1980 have provided a highly accurate set of tracking 
station coordinates. This independently derived solution has provided a 
means of independently verifying the Lageos derived values and has thus 
placed bounds on the overall precision relative to the center-of-mass of 
the coordinate system. 

Comparisons of the STARLETTE derived tracking station coordinates with 
those derived from LAGEOS at GSFC (GSFC SL5) and at the University of 
Texas (UT8081 A-M) are presented in Table 1. In these comparisons, 
translation and rotation parameters have been adjusted in order to removed 
any systematic differences in the coordinate system definitions. The 
AX and AY translation parameters were small (<10 cm) however an 
unexplained AZ shift of about 70 cm is noted between the LAGEOS 
solutions. The rms agreement of the Starlette derived values with the 
LAGEOS solutions is generally better than 20 cm. 

The recovered value for GM is 398600.43 ± .01 km^/sec^ which is in good 
accord with the value of 398600.436 recovered by Smith, et al. (1983) 
from 1980 LAGEOS data. This indicates the overall scale control in the 
solution is good and that the potential model is reasonably complete. 

From our previous work on STARLETTE presented in Marsh and Williamson 
(1° 7 8) and in Fel sentreger, et al . (1979), we were aware of the fact that 
STARLETTE is significantly perturbed by tides. The earth and ocean 
tidal effects on satellite position in a typical 5-da.y arc of STARLETTE 
are substantially over a meter. Given a sparse data set, there is little 
hope of attaining 10 cm station positions or baseline determinations 
without appropriately modeling this effect. At the time we designed this 
current solution, our software system at GSFC, GE0DYN, had recently been 
enhanced with the capability to model the tides as individual 
constituents. This modeling, developed by 0. Christodoulidis, pro- 
vides for linear responses at each specific frequency in the tide gene- 
rating potential and is expressed in terms of the osculating ecliptic 
elements of the disturbing body and the equatorial position of the 
desired point on the earth's surface. 


The approach we were able to use in this study is to recoyer tidal 
parameters directly in the simultaneous least squares data reduction 
along with the other geodetic parameters— tracking station locations, 
geopotential coefficients, and the like. This technique is possible only 
because of the improvement in the quality of the laser tracking data, in 
our knowledge of the geopotential, and in our otheV improved modeling 

related to centimeter satellite geodesy and geodynamics. 

' ^ 

Table 2 presents the STARLETTE recovered tides. These tides are the full 
tidal effect relative to the Wahr free space potential, and the Woolard 
nutation series coordinate system. Where significant sideband terms 
existed (typically those differing by the period of the lunar node from 
the dominant central term), the analytic theory of Cartwright and Taylor 
(1971) and the later Cartwright and Edden (1973) have been used to 
determine the appropriate value at the central frequency. The basic 
assumption in this projection is that the earth's response function 
(admittance) should be approximately the same at nearby frequencies. 
Otherwise, the values are as recovered. Note that the diurnal tides are 
appropriate to the coordinate system definitions described which incor- 
porate the Woolard series nutations. The uncertainties in these values 
range from about 0.5 cm in C 20 to about 0.2 cm in C22* 

SIGNIFICANCE 

Analyses of the 1980 Starlette tracking data have: 

0 independently confirmed the relative positioning of the LAGEOS 
recovered tracking station locations to an accuracy in each 
coordinate of 20 cm. 

0 illustrated a sub-meter bias in the AZ network center of 
origin in the satellite laser systems. 

0 demonstrated that the tidal response of the earth has been mapped 
to better than .01 in effective k 2 for the entire tidal frequency 
range from the direct analyses of the laser tracking data. 
This includes the determination of long period tides which are 
fundamental to UT1 variations. This mapping substantially 
agrees with a combination of the Schwiderski ocean tides and 
the Wahr body tide. 

0 demonstrated that the precision laser tracked satellites are 
quite sensitive to nutation error and in fact may be a good way 
to measure it. 

0 demonstrated that STARLETTE can be used to recover reasonable 
polar motion values after suitable corrections are made for 
nutation errors. 


FUTURE EMPHASIS 

Additional analyses using more recent laser tracking data recorded during 
1981 to 1983 are in progress. Techniques to better account for the effect 
of using the Woolard nutation series as versus using the modern Wahr 
series are also in progress. These analyses are expected to provide a 
further refinement of the present parameters and further insight into 
precision orbit determination at the decimeter level. 
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Table 1. Tracking Station Coordinate Comparison 
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Translation Parameters 
AX AY az 

X 

Rms of Fit 
Y 

Z 

Starlette vs. 
UT8081 A-M 

-9 cm. 

-4 cm. 

-43 cm. 

19 cm. 

25 cm. 

19 cm. 

Starlette vs. 
GSFC SL5 

8 

4 

29 

17 

on 

tuv 
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6SFC SL5 vs. 
UT8081 A-M 

2 

-8 

-72 

6 

10 
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7063 Greenbelt, Maryland 

7086 McDonald Observatory, Texas 

7090 Yarragadee, Australia 
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7114 Owens Valley, California 

7115 Goldstone, California 

7907 Arequipa, Peru 

7943 Orroral Valley, Australia 
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Table 2 

Tidal Coefficients from 1980 Starlette Data Analyses 
Relative to Wahr Earth Potential and Woolard Nutation 


Second Degree 


Fourth Degree 



Amplitude — C 2m 
(cm) 

Phase* - e 2m 
(deg) 

M 

065.455 

1.99 

42.7 

M f 

075.555 

1,70 

175.4 


085.455 

1.47 

149.6 

Q 1 

135.655 

2.14 

282.1 

Ot 

145.555 

3.32 

255.5 

M-j 

155.655 

.37 

284.6 


162.556 

1.14 

35.7 

Pi 

163.555 

.92 

106.1 

5<1 

165.555 

3.28 

32.2 

Jl 

175.455 

1.34 

75.9 

00 ! 

185.555 

4.56 

109.9 

2N 2 

235.755 

.05 

159.9 

n 2 

245.655 

.76 

232.9 

m 2 

255.555 

2.57 

229.2 

t 2 

272.556 

.27 

13.4 

s 2 

273.555 

.78 

302.4 

k 2 

275.555 

.83 

293.4 







Doodson (1921) convention used for Phase Definition 



















GEODYN PROGRAM SYSTEMS DEVELOPMENT 


B.H. Putney 


OBJECTIVES 

The purpose of the Geodyn Orbit Determination and Parameter Estimation, 
the SOLVE and ERODYN Programs is to recover geodetic and geophysical 
parameters from satellite and other data in a state-of-the-art manner. 

BACKGROUND 

In 1971 the NONAME and GEOSTAR programs were combined to create the GEODYN 
program. The SOLVE program was created at the same time. A few years later 
the ERODYN, error analysis program was written. The philosophy of the 
development of the software system has been maintenace of computer effi- 
cient, well -structured software, with appropriate orbit, earth and 
numerical models, using precise satellite measurement modeling and effi- 
cient numerical models, and performing careful benchmark procedures. 
This care has paid off in the production of several GEM'S (Goddard Earth 
Models), precision station locations, improved tidal, GM, and polar 
motion and earth rotation values and consistent baselines. Careful 
usage, analysis, and modeling using laser, altimeter and other 
satellite data from LAGEGS, SEASAT, STARLETTE, GEOS and BE-C satellites 
as well as many others has made these accomplishments possible. 

RECENT ACCOMPLISHMENTS AND SIGNIFICANCE 

As part of the Crustal Dynamics Project a significantly improved gravity 
field, pole positions, earth rotation, GM and baselines have been deter- 
mined. The new solid earth and ocean tidal model and observation normal 
point technique are being used and evaluated. 

In support of GRM, several mission successful simulations were performed 
using spherical harmonics, geoid heights and gravity anomalies. Numerical 
experiments with block interpolation,* more accurate block integra- 
tion, Fourier analysis of geopotential parameters and Householder trans- 
formation techniques for least squares estimation, have been studied. 
In preparing for the MERIT Campaign the DE 118 JPL sun, moon and plane- 
tary ephemerides are available in GEODYN 8210 version. This provides us 
with improved planetary positions and the Wahr nutation series. In 
addition, the ocean loading displacement is modeled in that same program. 

Conversion of GEODYN and SOLVE to the Control Data Corporation Cyber 205 
has been accomplished this year. Both programs are being optimized for 
that computer. GEODYN II is a rewrite of the GEODYN program. When 
completed it will execute on the Cyber 205 and on IBM compatible 
computers. Initial timings on the sequential GEODYN program show a 20 
percent improvement over the IBM 360/95 computer. Solve is running 
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three times faster on the Cyber than on the IBM 3081. We are pleased 
with the progress on the new computer. Modifications to make this a 
flexible tool, by allowing multiple satellites (GPS) to be handled in a 
single arc, processing a pass of observations simultaneously, allowing 
variation in numerical stepsize for high and low degree and order spheri- 
cal harmonic coefficients, making easier additions of new adjusted para- 
meters, restructuring of observation and input files to allow for more 
information and precision, are designed into the new program. The SOLVE 
program is currently being vectorized but its basic structure appears to 
be appropriate for the Cyber. However, modifications of Input /Output 
techniques are required for optimal use of the Cyber and these modif- 
ications are underway. Several new versions of GEODYN, SOLVE and ERODYN 
were released this year. Versions of the program were requested and 
delivered to England, Germany, Japan, the Netherlands, Australia, JPL 
and Ohio State University. 

FUTURE EMPHASIS 

The conversion and optimization of the GEODYN and SOLVE programs on the 
Cyber 205 will be completed. It is hoped to begin the conversion of the 
ERODYN program as well. At this time it is very important to thoroughly 
check out the GEODYN program models, i.e., albedo, relativity, plate 
motion, tidal model, coordinate system. Reexamining of approximating 
assumptions for correctness and for now required accuracies, of formu- 
lation and implementation Is now necessary. Additional changes for the 
MERIT campaign such as changing the reference year from 1950 to 2000 
need to be implemented. It is desirable to model and solve for a set of 
spherical harmonic coefficients for the sea surface topography and to 
incorporate a range difference data type to be used with laser data. 

The software will be modified to meet the needs of the Topex mission as 
they are defined and required. Additional GRM modifications will be 
made as defined and funded. 

Every effort will be made to continue to mold the software in a timely 
manner to support the needs of the scientific community determining 
earth models. 


I 
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CHAPTER 4 


SEA SURFACE TOPOGRAPHY 


OVERVIEW 


The overall goal of this area of research is to apply remote sensing 
data to studies of dynamic ocean and earth processes. The primary data 
types are satellite altimetry data, orbit perturbation data and gravity 
data. These data are analyzed in conjunction with surface observations, 
for example, ocean temperature and density measurements. 

Global sets of satellite altimeter data have been collected during two 
satellite missions: The GEOS-3 mission (1975 to 1978) and the Seasat 
mission (1978). 

Research programs are underway in the following areas to analyze and 
interpret the altimeter data: (I) the computation of highly precise 
regional mean sea surfaces which can be used for studying the oceanic 
mesoscale variability and to perform detailed analyses of the earth's 
crustal structure; (2) the study of lithospheric flexure seaward of 
subduction zones to obtain a better mechanical model of oceanic 
lithosphere; (3) to develop and test an interpolation technique which 
allows accurate extrapolation of tidal height fields in the ocean basins 
by making use of selected satellite altimetry and conventional tidal 
gauge measurements. 

Contributors to this chapter include J.G. Marsh, D.C. McAdoo and 
B.V. Sanchez. 
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REGIONAL MEAN' SEA SURFACES BASED UPON 
GEOS-3 AND SEASAT ALTIMETER DATA 


James G. Marsh 


OBJECTIVE 

The objective of this work has been the computation of highly precise 
regional mean sea surfaces based upon crossing arc adjustments of the 
combined GEOS-3 and SEASAT altimeter data sets. These surfaces have 
provided a basis for studying the oceanic mesoscale variability and for 
detailed analyses of the earth’s crustal structure. 

BACKGROUND 

The GEOS-3 (1975-1978) and SEASAT (1978) altimeter experiments provided 
very dense coverage in the N.W. Atlantic, the Gulf of Mexico and the 
Bering Sea. The number of available passes ranged from 558 in the Gulf 
of Mexico to 1284 in the N.W. Atlantic. The GEOS-3 data were smoothed 
over a frame to provide data at a rate of 1 point/second. Only intensive 
mode data were used, and gross data outliers were eliminated in the 
smoothing process. The SEASAT data used were the JPL distributed 
Geophysical Data Records (GDR) also at the rate of 1 point/second. 

The precision of the GEOS-3 altimeter measurements is about 25 cm and 
that of the SEASAT altimeter measurements is about 3-10 cm. However, 
in order to achieve precision of this order in the computed sea surfaces, 
other sources of error, chiefly the radial orbit error of the spacecraft, 
which may be on the order of a few meters, must be removed. This is 
done in the crossing arcs least squares adjustment process. 

The satellite radial orbit error is primarily due to the errors in the 
gravity field model. Several studies have shown that most of the radial 
orbit error has frequencies close to once per rev. Thus, over distances 
of thousands of kilometers, the radial orbit error can be well repre- 
sented by a linear function of time, and can therefore be removed by 
the least squares adjustment process if the parameters of a linear 
trend are solved for. Any other long wavelength sources of error, for 
example errors caused by mis-modeling of the ocean tides, will also be 
removed in the adjustment process. The adjustment process has been 
extensively described elsewhere (Marsh, et al., 1980). 

The large amount of data in these areas coupled with the increased 
constraints provided by the combination of data from two orbital 
inclinations has provided the basis for a detailed precision adjustment 
of the data. The precision of these surfaces is estimated to be on the 
order of 15 cm with a horizontal resolution of 25 km. 
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RECENT ACCOMPLISHMENTS 


Figure 1 shows the sea surface contour map of the Northwest Atlantic 
derived from the SEASAT and GEOS-3 altimeter data. Several features of 
the ocean floor topography are clearly reflected in the sea surface 
topography. Bermuda is clearly visible, and the Ibng minimum in the 
center corresponds to the Hatteras abyssal plain. To the west, the 
Blake Bahama Outer Ridge and the Blake Escarpment are visible. 

Figure 2 shows the sea surface computed in the Bering Sea using SEASAT 
and GEOS-3 altimeter data. The sur f ace clearly shows where the shallow 
ocean over the continental shelf (in the northeastern portion of the 
map) is divided from the deeper ocean over the Bering Abyssal Plain in 
the southwest of the sea. Also clearly visible are the ocean bottom 
ridge at about 170° longitude, and the feature known as Bowers Bank 
(the ci rcular feature below latitude 55° and between 175° and 180° 
longitude). 

The details of this work are presented in a report by Marsh et al . 
(1983). 

SIGNIFICANCE 

Except for mean dynamic topography, the altimstric surfaces presented 
in this paper are equivalent to the geoid. It should thus be possible 
to use these surfaces to remove the geoid signal from synoptic altimeter 
data in order to detect transient mesoscale features. Of course, one 
must proceed with caution for a quasi -permanent current such as the 
Gulf Stream, since both synoptic and time-averaged altimeter data 
contain a Gulf Stream signal. In contrast, gravimetric reference 
surfaces can be used to detect permanent circulation features. This is 
demonstrated in Figure 3 where we have differenced the Atlantic 3-year 
mean surface with the detailed 5' geoid model of Marsh and Chang (1978). 
The Gulf Stream can be clearly seen as a 1-m step running eastward from 
Cape Hatteras. South of the Gulf Stream is a broad area in which the 
two reference surfaces show agreement to within approximately 20 cm. 
Even though this area is populated by cold rings with sea height 
signatures of 50-100 cm, the transient nature of these features causes 
them to be considerably smoothed in the 3-year mean surface. Large 
differences occur in the northeast due to the New England seamounts and 
in the southwest due to the Bahamas. These differences are attributable 
principally to errors in the detailed geoid model. 

FUTURE EMPHASIS 

Future work in this area will be oriented toward further definition of 
mesoscale variability from crossover difference statistics, the detec- 
tion and monitoring of individual cold and warm core rings using the 
mean sea surfaces as a reference and toward the orientation of these 
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highly precise mean sea surfaces into an absolute center-of-mass coor- 
dinate reference system. 
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Figure 1. Mean Sea Surface Based Upon SEASAT and GEOS-3 Altimeter Data, 1975 - 1978 
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SEASAT OBSERVATIONS OF LITHOSPHERIC FLEXURE 
David C. McAdoo 


OBJECTIVE 

The individual lithospheric plates move about the Earth's surface as units 
which remain relatively rigid on geological time scales. It is of con- 
siderable importance for geophysicists to know as much as possible about 
the actual strength of the lithosphere including the dependence of 
strength on age. This study of lithospheric flexure seaward of subduc- 
tion zones seeks a better mechanical model of oceanic lithosphere. The 
study is based primarily on SEASAT altimeter observations. An improved 
understanding of forces acting at subduction zones is also sought. 

BACKGROUND 

The Outer Rise seaward of deep ocean trenches is generally attributed to 
flexure of the oceanic lithosphere prior to subduction. Various models have 
been proposed to describe this flexure. They all include a mechanically 
strong lithosphere — loaded from the arcward side of the trench axis — overlying 
a fluid-like asthenosphere. Simplest among them is the uniform elastic plate 
model which was first applied to Outer Rises by Hanks (1971), and subsequently 
by a number of others. 

More complex models which include a variable rheology with depth have also 
been developed. These complex models produce more realistic estimates of 
stress but, aside from this, the elastic model hs been shown (Bodine, et al., 
1981) to be adequate for describing flexure of the oceanic lithosphere. 

Quite naturally these lithospheric models have been tested almost exclusively 
on bathymetric profiles. However, the Outer Rise also has a clear gravita- 
tional expression. Watts and Talwani (1974) showed that gravity anomalies 
over the Outer Rijse tend to mimic— in form— the bathymetry. It should, 
therefore, not be surprising that a geoidal counterpart of the Outer Rise 
exists and is evident in the SEASAT altimeter observations of the marine 
geoid. This study extends the popular elastic plate model to include a 
predicted geoid signal and then tests this model directly on SEASAT observa- 
tions of the marine geoid. One anticipated result from this study is a 
clearer understanding of the dependence of the lithosphere's effective elastic 
thickness (EET) on age. 
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Using bathymetric data Bodine, et al . (1981) concluded that— in the absence 
of horizontal loading— the effective thickness, T e , increases with age accor- 
ding to « C * age 1 '^. Their results suggest the thickness (EET) is, 
in several regions, significantly less than the value predicted by this 
thickness-age relation; they attributed this discrepancy to horizontal com- 
pres'sive loads acting across these trenches. 

RECENT ACCOMPLISHMENTS 

A paper, McAdoo and Martin (1983) has been completed. In this paper, segments 
of 102 SEASAT altimeter passes over six trench /Outer Rise complexes have been 
analyzed. For each trench the number of passes analyzed is a follows: Middle 
America, 14; Philippine, 12; Aleutian, 23; Mariana, 12; Izu-Bonin, 10; Kuril, 
30. Thes sea surface height profiles obtained from these passes will be 
treated as geoid height profiles. 

After long-wavelength component of the geoid is removed, these observed geoid 
height profiles— i .e., the SEASAT passes are then compared to geoid profiles 
predicted by the elastic model. As an example, one pass over Mariana Trench 
is shown in Figure 1. 

Using the method of least squares, model parameters are adjusted and estimated 
for each SEASAT pass. For each pass an estimate of flexural wavelength is 
obtained. A mean flexural wavelength is then obtained from the ensemble of 
SEASAT passes over each trench. From this mean flexural wavelength, an 
effective elastic thickness can be derived. The six crosses in Figure 2 
represent the effective elastic thickness versus age of the subducting oceanic 
lithosphere for each of the six trenches studied. The solid curves represent 
the thickness age law predicted by Bodine, et al . (1981) from thermal and 
experimental constraints. The agreement between this law and effective 
elastic thicknesses derived from SEASAT data is good; in fact it is better 
than the agreement obtained in previous studies (see Bodine, et al., 1981) 
which use bathymetric data— not SEASAT data— to estimate thicknesses (EET). 
The lack of agreement obtained in these previous studies was generally attri- 
buted to unmodeled horizontal compressive forces acting across certain 
trenches. This work is described in McAdoo and Martin (1983); see also 
McAdoo, et al . (1982). 
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SIGNIFICANCE 


SEASAT altimeter data provide important new information on the deformation 
of the lithosphere seaward of deep ocean trenches. These data, when com- 
pared with a mechanical model, produce estimates of the effective elastic 
thickness of the lithosphere which are more consistent from profile to 
profile in a given region than those derived from bathymetry. For the six 
trenches which we studied, the effective elastic thickness derived from 
altimetry increases with the square root of lithospheric age as predicted 
by Bodine, et al . (1981). Significant horizontal compression across 

these trenches is not indicated. 

FUTURE EMPHASIS 

The potential utility of SEASAT and other satellite altimeter data for 
studying dynamics of the lithosphere and sub-lithosphere is only beginning 
to be exploited. This particular study of outer rises will be extended to 
include all the deep ocean trenches. Mechanical structure of the litho- 
sphere can be similarly studied in mid-plate regions. These types of 
studies which refine our models of the lithosphere are necessary if we are 
to isolate, with confidence, dynamical effects arising beneath the lithos- 
phere such as small scale convection in the mantle. 
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Figure 1. A residual SEASAT height profile compared with best fitting 
result from the elastic plate model. 
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Figure 2. Effective elastic thickness (EET) derived from SEASAT data (crosses) 
over six trenches (Middle America (MD), Aleutian (A), Philippine (P), Kuril (K), 
Mariana (M) and Bonin (B). EET's are plotted versus age. Bodine, et al.'s (1981) 
age-thickness relation is shown as curves. 







AN OBJECTIVE ANALYSIS TECHNIQUE FOR EXTRAPOLATING TIDAL FIELDS 

Braulio V, Sanchez 


OBJECTIVES 


1. ' To develop and test an interpolation technique which allows 

accurate extrapolation of tidal height fields in the ocean 
basins by making use of selected satellite altimetry measure- 
ments and/or conventional gauge measurements. 

2, To develop a dynamical model of the ocean tides. 


3. To test the technique initially for the M2 tide in a small 
water body such as Lake Superior by theoretically computing 
the tide and to compare the results with available tidal 
data. 

4. To test the interpolation technique in Lake Superior in 
order to ascertain its capabilities for further application 
in the ocean basins. 

BACKGROUND 

The theoretical foundation of the method used in the calculations is 
based on Proudman's theory (1918) as reformulated by D.B. Rao (1966). 
The theory provides a formulation for calculation of gravitational 
normal modes and rotational normal modes of irregularly shaped basins 
with realistic bathymetry. The method has been extended to compute 
the forced solution (M2 tide) in Lake Superior to serve as a test to 
evaluate the objective analysis procedure. The objective analysis 
uses orthogonal functions that formed the basis to represent the water 
level fluctuation field in the free and forced calculations. The 
available tidal amplitudes and phases are expressed as a linear 
combination of these functions. Since data points are usually avail- 
able only at irregularly spaced points* instead of using the orthogo- 
nality property of these functions to determine the expansion coeffi- 
cients, they are determined through a least squares procedure. 
In the computation of the forced solution we include the effects of 
the yielding of the solid earth to tide generating forces. 
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RECENT ACCOMPLISHMENTS 


a. The M2 tide solution 

Figure 1 shows the amplitude and phase fields for the M2 tide in 
Lake Superior. The maximum amplitudes* which are close to 20 mm, 
occur at the western end of the basin. The phase values obtained 
by Mortimer and Fee (1976) from analysis of records of water levels 
are given in parentheses. The amplitudes values were provided 
by C.H. Mortimer (personal communication) from the original 
calculations on the water level data. Also given are the phase 
and amplitude values obtained in this investigation. The latter 
are the values obtained at the grid points closest to the stations 
in Mortimer and Fee’s analysis. The contours of equal amplitude 
and phase are given by the solid and dotted lines respectively. 
The solution was obtained by including the first SO eigenfunctions 
from the velocity potential and stream function solutions. The 
forced solution then contains 80 coefficients, the magnitudes of 
the coefficients indicate the relative importance of the different 
characteristic functions in accounting for the variance of the 
total solution. Most of the power is contained in the first coef- 
ficient. The results for the first 12 coefficients are given in 
Table 1. An examination of Figure 1 shows excellent agreement 
between the theoretically computed tidal response and observed 
amplitudes and phases. The theoretical simulation also serves as 
a guide in choosing the most energetic modes to be used in the 
objective analysis at any order of truncation of the repre- 
sentation. 

b. Interpolation Technique Evaluation 

The data locations selected for the analysis are the ten water 
level stations analyzed by Mortimer and Fee. The tidal ampli- 
tudes and phases at these ten locations were extracted from the 
theoretical simulation. An objective extrapolation of the tidal 
field is carried out with these 10 data points using 2,4, 5, 8, 
and 10 most energetic modes listed in Table 1. A comparison of 
the objectively analyzed amplitudes and phases with the "true" 
values at all grid points is shown in Table 2. The table lists 
the percentage of total grid points where the analyzed and true 
values of amplitude and phase agree with each other within the 
limits indicated and for the number modes chosen in the objective 
analysis. There are a total of 227 points within the basin. 
Also given in the table are the rms errors of the amplitude over 
the entire basin and only at the data points, as well as the 
ratio of the potential energy of the analyzed and true fields. 
It is seen from Table 2 that the method of analysis yields a 
tidal height field that is within 3% of the "true" potential 
energy with just 2 basis functions. As the number of functions 
are increased to 5, the energy difference is within 5%. For the 
ten selected data points, the fit with 5 functions appears to be 
the best. The tidal pattern using 5 coefficients is shown in 
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Fig. 2.' Comparison with Fig. 1, the theoretical solution obtained with 80 
coefficients, shows excellent agreement between the two fields. As the 
expansion coefficients are increased to 10, which represents an exact deter- 
mination of the values of the coefficients, the overall fit at the tidal 
pattern deteriorates. Obviously, the maximum number of coefficients that 
would result in a best fit over the entire basin depends on the number and 
distribution (in space) of the available data points. The combination of 
theoretical simulation and objective analysis also helps in choosing loca- 
tions for additional observations that would contribute most to the res- 
olution of the tidal field in an ocean basin. 

SIGNIFICANCE AND FUTURE EMPHASIS 

A method has been described to objectively analyze available tidal ampli- 
tudes and phases to map the tidal constituent over an entire basin. The 
method consists of constructing a set of orthogonal functions that are 
determined by the geometry and topography of the basin. These functions 
are linearly combined and the expansion coefficients are determined from 
available data in such a way as to minimize the mean-squared error by a 
simple matrix inversion. Since the orthogonal functions used are charac- 
teristic of the given basin, they would reflect the basin properties in a 
sense and hence would be capable of giving the best fit for the overall 
tidal field in the basin. Using Lake Superior as a test case, a theore- 
tical simulation was done for the M2 tidal amplitudes and phases in the 
basin. This theoretically simulated field has been used as a "true" solu- 
tion to evaluate the recovery of the total tidal field in the basin from 
application of the analysis procedure to data taken at selected points. 
The agreement between the analyzed and the "true" fields is found to be 
excellent. The theoretically calculated M2 tidal heights and phases also 
agree very well with those obtained from analysis of water level data by 
Mortimer and Fee. The method is now being extended to ocean basins. The 
results of the Lake Superior investigation have been submitted for publi- 
cation by Sanchez, Rao and Wolf son. 
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Figure 1. Theoretical solution for the M2 tide in Lake Superior (80 coefficients). Contours of equal 
amplitude (mm) and equal phase (degrees) are given by the solid and dotted lines respectively. 
The values from Mortimer and Fee's data reduction are given in parentheses. 
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Figure 2. Solution for the M2 tide based on objective analysis using 5 coefficients corresponding to the 
five most powerful modes. Contours of equal amplitude (mm) and equal phase (degrees) are given by 
the solid and dotted lines respectively. The values from the theoretical solution (80 coefficients) 

used as data are given <n oarentheses. 
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Table 1 

Power Spectrum for the M2 Tide 


Mode 

Power (%) 

Period (hours) 

1 

94.62 

7.83 

2 

1.74 

4.42 

3 

0,70 

3.75 

4 

0.95 

3.16 

5 

0.96 

2.96 

6 

0.05 

2.57 

7 

0.03 

2.40 

8 

0.35 

2.32 

9 

0.05 

2.20 

10 

0.15 

1.91 

11 

0.09 

1.82 

12 

0.10 

1.65 


Coeffi- 

cients 


Table 2 

Comparison between Objectively Analyzed and Theoretical Fields 


Amplitude 



<5% <10% <20% 


23.8% 
26.9% 
42.3% 
25.1% 
% 



<20% 

O 

in 

V 

<10° 

| 

<20° 


RMS 

Error 

(mm) 

P T 

64.8% 

46.7% 

68,7% 

83.3% 

1.0295 

1,2006 

76.7% 

64.8% 

85.0% 

96.0% 

1.0749 

1.2059 

87.2% 

61.2% 

86.8% 

93.4% 

0.9521 

0,8838 

55.5% 

59.9% 

87.2% 

96.5% 

1.0872 

3.1046 

36.1% 

24.7% 

39.6% 

58,6% 

3.9678 

12.185 




























CHAPTER 5 


ADVANCED STUDIES 


OVERVIEW 


The advanced studies activities considered by the Geodynamics Branch provide 
a sound basis for the planning, design, and develop^^t of sensor and space 
mission concepts. Studies have been conducted based on extensive current 
research experience and simulations. The utility of the data expected from 
these systems is evaluated in terms of predicted contributions to program 
objectives and user applications. 

Currently, the horizontal resolution of satellite-determined gravity field 
models is of the order of 500 km. This level of resolution must be increased 
to meet requirements for solid-earth and ocean studies. A gravity field 
mapping mission, namely The Geopotential Research Mission (GRM) is planned 
for the latter part of this decade to map the eartfrVs gravity and magnetic 
fields. GRM consists of two spcecraft in the same circular/polar orbit 
with a separation distance varying from 1° to 6° and a nominal orbit altitude 
of 160 km. From Doppler measurements made between the two satellites with 
a precision of 1 micrometer per second it is the goal for GRM to map the 
earth's gravity field to an accuracy of 1 mil li gal and the geoid to an 
accuracy of 5 cm with a horizontal resolution of 100 km. 

Two GRM simulation studies are presented in this chapter. The studies 
describe (1) the gravity field signal characteristics of the GRM Doppler 
sensor and (2) the computer simulations which were performed to assess the 
capability of GRM to map the gravity field of the earth for geophysical and 
oceanographic applictions. Another new sensor is discussed in this chapter 
which provides a method of precise isitioning of a large number of points 
on the earth's surface in a short period of time. That is, a measure of 
the relative location of geodetic markers from an aircraft. Such a system, 
"The Airborne Laser Ranging System (ALRS) is a potential sensor for per- 
forming rapid (few days) mapping of local crustal motion in seimsically 
active areas. 

In a somewhat different vein, an analysis of soil moisture based on measure- 
ments of the earth surface temperature has also been under study. Although 
this study is most applicable to problems in non-renewable resources, the 
coupling of soil moisture and atmospheric conditions to geodynamic processes 
is also under investigation. 

Contributors to this section are: T.L. Felsentreger, F.J. Lerch, W.D. Kahn, 
and J.E. Welker. 
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GRAVITY FIELD SIGNAL ANALYSIS STUDY FOR A GEOPOTENTIAL RESEARCH MISSION 

(GRM) 


OBJECTIVE 


Theodore L. Felsentreger 


The objective of this report is to discuss the gravity field signal 
characteristics of the SST Doppler sensor to be flown on GRM. 

BACKGROUND AND RECENT ACCOMPLISHMENTS 


The Geopotential Research Mission (GRM), formerly known as GRAVSAT/MAGSAT, 
has been proposed for the mapping of the Earth's gravity field. It is a 
candidate for a NASA new start in the near future. GRM is to be composed 
of a pair of surface-force compensated satellites in identical polar 
orbits, spaced approximately 3° apart at a nominal altitude of 160 km 
with a system precision of 1 ym/sec. 

The perturbation by the gravity field induced on the intersatellite 
Doppler signal is given by the following relationship: 


where 


5 ( r i) 
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r d\p 


dt 


(1) 


T = disturbing potential 

a = Earth's mean equatorial radius 

h = satellite orbit altitude 

r s a+h 

t|) s true anomaly 

St= time separation of two satellites in identical orbits. 


The form of the disturbing potential most useful for formulating (1) is 
given by the formula 
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Ag = mean gravity anomaly 
da = element of area. 
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Using (1) and (2), velocity profiles resulting from the perturbation 
produced by a (1° xl°) 1 mgal gravity anomaly as sensed by the two 
satellites are shown in Figure 1— results are shown for 160 km and 200 km 
satellite altitudes. The amplitude of the gravity signal drops off by a 
factor of 1.5 when the altitude increases from 160 km to 200 km. By 
extrapolation the amplitude is further decreased by * about a factor of 3 
when the orbital altitude is increased to 250 km. Thus, the amplitude of 
the measurement drops off as the altitude increases, to the point where 
the measurement amplitude is insignificant at the 1 mgal level for alti- 
tudes above 200 km. By way of contrast. Figure 2 shows the SST signal 
generated by GRM for one full orbit, under the influence of a full 
(36,36) spherical harmonic gravity field (GEM 10B). 

The ability to resolve the signal from a (1° x 1°) 1 mgal gravity anomaly 
for two spacecraft whose orbits are displaced from the anomaly by 1° and 
5° is shown in Figure 3. For a 1° displacement, better than 60% of the 
gravity signal is sensed; however, for a 5° displacement, only 10% of the 
gravity signal is sensed. This implies that the gravity anomaly signals 
from adjacent blocks drop off drastically with increased distance from 
the satellite ground track. Nevertheless, it is significant that remote 
gravity anomalies do contribute appreciably to the overall gravity signal. 

FUTURE EMPHASIS 

Studies of this type will be continued, along with rapid error analysis 
and full-scale simulation studies. 
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Figure 1. GRM Anomaly Sensing at Two Altitudes. 
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Figure 2. GRM Full (36, 36} Gravity Field Sensing 
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GRAVSAT SIMULATIONS FOR GEOPOTENTIAL ANALYSIS 
Francis J* Lerch 

OBJECTIVES 

The objective of this work is to assess through computer simulations the capa- 
bility of the Geopotential Research Mission (formerly Gravsat for improving 
the gravity field of the earth for geophysical and oceanographic applications. 
Objectives of the mission are 3 mgal for gravity anomaly for geophysical analysis 
and 10 cm for the geoid for ocean circulation with resolution down to 1° half- 
wavelength. 

BACKGROUND 

The Geoptential Research Mission (GRM) through use of a dedicated satellite 
will map the gravity field fine structure over the earth and thus provide 
knowledge and understanding of (1) the origin and structure of geological 
features on the earth's surface; (2) the mechanical properties of the earth's 
lithospheric plates and the forces which drive their motion; and (3) the large- 
scale circulation of the oceans and major current systems by the determination 
of an improved ocean geoid. 

The mission configuration for the GRM spacecraft consisted of two low-altitude 
satellites at 160 km in a drag free system separated by 300 km (used for 
simulations) with satellite-to-satellite doppler tracking (±1 micron/sec) 
surveying the globe in polar orbits. Error analysis studies have been made for 
GRM (see references). However, the main benefit of the present technique is 
that a simulation is performed from numerical integration of the force equations 
where a set of exact GRM observations are computed from an assumed gravity 
model (ground truth) which is to be recovered from the study. This method is 
free of a number of different approximations made in other studies but is 
generally costly in computer time and hence must be limited in scope. 

RECENT ACCOMPLISHMENTS 

Computer simulation results for the Geopotential Research Mission (formerly 
Gravsat) have been obtained which support accuracy objectives of the mission, 
namely 3 mgals for gravity anomaly and 10 cm for geoid height with a half-wave- 
length resolution of 1 degree. The computer simulations were performed with 
satellite-to-satellite doppler tracking (+ 1 micron/sec) between two low alti- 
tude satellites at 160 km in polar orbits and separated by 300 km in a drag- 
free environment for a 6-month mission lifetime. Results from studies have 
been obtained as follows: (1) the recovery of spherical harmonics with poten- 

tial coefficients complete through degree and order 36 gave an accuracy of 2 
to 3 orders of magnitude better than current knowledge; (2) recovery of geoid 
height and gravity anomaly in a local region (50° LAT x 20° LONG) in the pre- 
sence of unmodeled errors external to the region gave an accuracy of 7 cm for 
geoid height and 4 mgal for the gravity anomaly for 1° x 1° blocks within the 
interior of the region; (3) the recovery of 1/2° x 1/2° blocks of geoid height 
and gravity anomaly ;in a local area containing a trench region gave, with use 
of eigenvalue analysis, commission errors of 5 cm for geoid height and 3 mgal 
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for gravity anomaly for 1° x 1° block sizes; and (4) a Fourier analysis over the 
trench region in (3) confirmed the results of this study which satisfy accuracy 
goals of GRM. Table 1 gives a further summary of these results. 

SIGNIFICANCE AND FUTURE EMPHASIS 

Results obtained so far are encouraging in that they satisfy mission goals but 
they are not complete and need to be extended. Future analysis should be 
conducted to extend the spherical harmonics to high degree (70 to 100) to reduce 
the aliasing effect on the solution due to the truncation of the harmonics. 
Because of computer time limitations it is necessary to truncate the spherical 
harmonic series. The residual potential including the smaller wavelengths of 
the truncated harmonics will then be recovered from localized functions by 
combining steps (2) and (3) above to cover a considerably broader region than 
the 50° x 20° area used in the present analysis. Also, because of computer time 
efficiency, the technique of deriving the potential first at satellite altitude 
and then downward continuing the result to the surface of the earth is being 
investi gated. 

TABLE 1. SUMMARY OF GRAVSAT SST SIMULATIONS 

1. Recovery of Spherical Harmonics to 36 x 36 

o 20 days of SST data only employed 
o 2 to 3 orders of magnitude improvement over GEM 9 
o Iteration will give another order of magnitude improvement 

o Aliasing error is significant so a high degree solution is required 

o Orbit parameters modeled with a priori constraints of 30 m in position 
3 cm/sec in velocity 

o Combined SST+NWL dopper tracking gave similar solution as above. 

2. Geoid Solution in a Local Region 

o Local solution of 1° x 1° geoid heights in a 20° x 50° region gave 
7 cm accuracy for the interior region 
o Error mostly due to aliasing effects of unmodeled residual potential 
exterior to the region 

o A buffer zone of a priori constraints between the interior and exterior 
regions was necessary to control the aliasing error 

3. Short -Wavelength Resolution in Trench Af^ii (Tonga Tr.) 

o 1/2° x 1/2° geoid heights and gravity anomalies simulated over a 
5° x 5° trench area 

o The recovery gave errors for 1° x 1° blocks of 5.0 cm for geoid 

height and 3.0 mgal for gravity anomaly by eigenvalue analysis. 
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THE AIRBORNE LASER RANGING SYSTEM 
Werner D. Kahn 


OBJECTIVES 

The Airborne Laser Ranging System (ALRS) is a potential sensor for performing 
rapid (few days) mapping of local crustal motion in seismically active areas. 
Through the use of covariance error analysis techniques, studies are performed 
to determine measurement precision, measurement rate accuracy, and geometry 
requirements for the ALRS to detect and monitor small earth surface motions. 

BACKGROUND 

The Airborne Laser Ranging Sytem is a proposed multi beam short pulse laser 
ranging system on board an aircraft. It simultaneously measures the distances 
betweeen the aircraft and six laser retroref lectors (targets) deployed on the 
earth's surface. Depending on the host aircraft and terrain characteristics, 
the system can interrogate hundreds of targets distributed over an area as 
large as 6 x 10 4 sq. kilometers in a matter of hours. Potentially, a total of 
1.3 million individual range measurements can be made in a six hour flight. 
The precision of these range measurements is approximately ±1 cm. These 
measurements are then used in a procedure which is basically an extension of 
trilateral on techniques to derive the intersite vector between the laser ground 
targets. By repeating the estimation of the intersite vector, strain and strain 
rate errors can be estimated. These quantities are essential for crustal 
dynamics studies which include determination and monitoring of regional strain 
in the vicinity of active fault zones, land subsidence, and edifice building 
preceding volcanic eruptions. 

RECENT ACCOMPLISHMENTS 

Simulations have shown that range measurements must be taken at two widely 
separated altitudes in order to strengthen the geometry sufficiently to recover 
baselines at the centimeter level in agreement with earlier results (1). Thus, 
in a typical mission, the aircraft (i.e., NASA's NP3A Orion) approaches the 
target grid at an altitude of 3.9 km as in Figure 1. After acquiring the first 
few targets, the instrument shifts to the tracking mode for the remainder of 
the mission. After overflying the rows of targets at 3.9 km, the aircraft 
climbs to its maximum cruise altitude (say 6 km) for a second set of passes 
over the target grid. 

The spacing between targets will depend on several factors including the 
scientific objectives of the mission, the aircaft altitude and the terrain 
limitations. For the NP3A aircraft, the spacing is nominally taken to be 7 km. 
At typical cruise velocities (i.e., 200 knots), laser range data to a given 
target is taken for approximately 130 seconds before the pointing system is 
commanded to acquire a new target. For a laser repetition rate of 10 pps, this 
corresponds to 1300 range measurements per target. Since a given target is 
common to a number of six target sets, and data are collected at two altitudes, 
approximately 7B00 measurements are typically made to each target. 
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Figure 2 shows the important of refraction errors and also illustrates that an 
extensive network of meteorological sensors is not required. Atmospheric 
measurements made at a single site collocated with a laser target within the 
ALRS grid, will significantly reduce the effects of atmospheric refraction upon 
ALRS baseline precision. In the estimation process, if the information from 
that single "met sensor" is used and the atmospheric parameters for the other 
ALRS targets are made part of the set of estimated parameters, a factor of about 
7 improvement in baseline precision is achieved. The figure also shows that 
the inclusion of additional meteorological sensors within the ALRS target area 
does not signficantly improve baseline precision. In the simulation it was 
assumed that the surface measurements of pressure and temperature were accurate 
to 1 mbar and 1.4°C respectively. The vertical variation in pressure was 
assumed to be determined by the hydrostatic equation. 

In Figure 3 the evolution of baseline precision is given for a series of randomly 
deployed laser targets. These targets are distributed (as shown in the inset) 
in a potential ALRS flight test region in the vicinity of Shenandoah, VA. The 
targets were located at approximately 17 first order survey monuments currently 
maintained by the U.S. Geological Survey. The simulation indicates that the 
random pattern does not significantly affect baseline precision as can be seen 
by comparing the results with those in Figure 2. The baseline error typically 
grows at the rate of 0.05 cm/km. For the simulation one meteorological sensor 
was located in the middle of the nap and the meteorological parameter at other 
sites were determined in the estimation process. 

The ALRS System has the potential capability to measure the relative horizontal 
and vertical locations of targets spaced tens to hundreds of kilometers apart 
to precisions of a few centimeters, it is thus an attractive candidate for use 
in a number of fundamental and applied crustal geodynamics measurements. Among 
these measurements are relative motions of tectonic plates across the plate 
boundaries and crustal determinations prior to and after major earthquakes. 

Baseline rate accuracies can be obtained from the baseline precision by using 
the covariance matrix from a least squares fit to a linear regression line. 
S.C. Cohen in (2) analyzed the level of precision with which crustal strain 
rates can be determined. The method used for the results given here follow his 
analysis. By assuming that the baselines are surveyed during a one-day period 
every 90 days, then as shown in Figure 4, baseline rate accuracies at the end 
of a ten-year observation period under the assumption of constant baseline 
motion can be determined for 5 km baselines to better than 0.02 cm/yr and for 
100 km baselines to about 0.1 cm/yr. By repeated measurements of the grid area 
within the seismically active region to the baseline rate accuracy level 
indicated, the ALRS could observe the deformation of plates at their boundaries, 
as well as general plate motion and attempt to "capture" a moderate sized 
earthquake, i.e., one of approximately magnitude 7 through repeated and densi- 
fied monitoring of the grid area once changes in the baseline rate are 
observed. 


FUTURE EMPHASIS 


Rigorous covariance error analysis studies for ALRS have been completed. 
Approval is sought for an ALRS demonstration experiment. Upon approval, com- 
puter software must be developed in order to perform simulations which are to 
verify the covariance analysis performed to date, also to test the closed 
loop data processing capability and to develop the experiment plan encompassing 
the operational aspects of the ALRS demonstration test. 
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ASSUMPTIONS : 

• ACFT ALTITUDES 
h MAX = 6 km 
h M I N = 3.9 km 

• ACFT SPEED: 0.1029 km/s (200 KNOTS) 

• NO. OF TARGETS: 17 (RANDOMLY DISTRIBUTED) 

• ALAS SYSTEM 
NOISE: ± 1 cm 
BIAS: t 1 cm 

• METEOROLOGICAL SENSORS 
a PRESS = ± 1 MBAR 

• DATA RATE 

10 MEAS/S ( - 127,000 ME AS) 


BASELINE DISTANCE (km) 


Figure 3. Airborne Laser Ranging System Baseline Precision for Random Laser Target Distribution. 
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SOIL MOISTURE FROM TEMPERATURE MEASUREMENTS AT THE EARTH'S SURFACE 


Jean E. Welker 


OBJECTIVE 

To determine soil moisture budgets at the earth's surface from soil and 
atmospheric temperature variations. 

BACKGROUND 

The soil surface and atmospheric interface can be considered a boundary 
between two dynamic systems. On the one hand, the atmospheric strata 
immediately above the earth's surface can be viewed as a number of layers, 
each responding individually and collectively to air flows and currents, 
to irradiation both from incoming solar radiation as well as reflected 
and emitted radiations from the ground surface below, and to rapid changes 
of moisture content in its various states and accompanying heat transfers. 
The soil immediately below the earth's surface can also be treated as layers 
undergoing heat transfers, temperature changes, and moisture flows. Both 
sides of the earth's interface, the 1 to 2 meter layer at either atmosphere 
above or soil below, are coupled in a complex manner to this surface interface 
by a variety of effects. 

From a one-dimensional view of temperature alone, variations at the earth's 
surface manifest themselves in two cyclic patterns of diurnal and annual 
periods, due principally to the effects of diurnal and seasonal changes in 
solar heating as well as gains and losses of available moisture. Beside 
these two well known cyclic patterns, a third cycle has been identified which 
occurs over a mesoscale period of roughly 3 to 10 days. This cycle, which 
occurs in diurnal maxima and minima extrema, is the subject of this current 
research. 

RECENT ACCOMPLISHMENTS 

A number of data sets were plotted and statistically analyzed in order to 
accentuate the existence and the characteristics of this mesoscale soil 
temperature extrema variations and its relations to other parameters. The 
cycle is initiated by a preciptation event, which bathes the soil and causes 
an immediate drop in the temperature extrema. This is followed by a tempera- 
ture recovery period of roughly a week in the absence of any more precipi- 
tation. In a 5-month, March through July period in northern mid-latitudes 
with no elongated periods of drought, the cycle can be expected to occur 
approximately one to three times a month. This temperature loss and re- 
covery can be seen in Figure 1, in which drops in maximum soil temperatures 
at depths of 5, 10 and 20 cm for the Tifton Meteorological Station in the 
State of Georgia for the period of March through July of 1979 are plotted. 
Other data plotted with maximum diurnal soil temperatures are maximum air 
and evaporation pan (4‘ in diameter) temperatures, recorded precipitation 
and pan evaporation. 
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This study was initiated using data from the State of Georgia, because of 
the contrasting severe drought and non-drought years which occurred in the 
late seventies-early eighties in that region. In contrast to the non- 
drought data in Figure 1, Figure 2 contains 1980 data which is characterized 
by a severe 27-day drought period in May and June with only a trace of rainfall 
on 1 day in this period. ' To eliminate soils as a critical factor in this 
study, Iowa station data with completely different soil characteristics were 
used in comparison with Georgian data. In order to broaden the contrast 
further, Soviet data was used in a good year/bad year comparison for the years 
1971 and 1972, respectively. Other earlier data sets form Eastern Europe were 
also included. 

All data sets are being analyzed in order to determine regional characteri- 
zation of drought conditions. A number of broad drought, descriptions were 
considered additionally in conjunction with the station temperature extrema 
data described above. 

Analyses of the data sets have been directed toward the behavior of tempera- 
ture extrema before and immediately after a precipitation event. With or 
without a drought condition, the diurnal minimum soil temperatures at 10 cm 
depth are highly correlated for 153 values, March through July diurnal values, 
with both the diurnal minimum air and evaporation pan temperatures. This is 
not to say that the 10 cm minimum soil temperatures and minimum air and evapo- 
ration pan temperatures have the same values without some offset, but only 
that they are highly correlated statistically, with a value greater than 0.85, 
for those cases tested for both drought and non-drought conditions. The 
presence of drought is manifested in a high correlation of diurnal maximum 
values rather than diurnal minimum values, for the data sets of the 10 cm soil 
temperatures, air temperatures, and evaporation pan temperatures. The 10 cm 
maximum soil temperatures and air temperatures move together under drought 
conditions because solar radiation at the ground during the hottest part of 
the day heats the soil and air together, with no energy loss to the evap- 
oration of soil moisture. After a precipitation event, in the period when 
the soil is hindered in regaining its original temperature because of heat 
absorption by the evaporation process, it is again the diurnal maximum rather 
than minimum soil temperatures that appear most sensitive. For all maximum 
and minimum temperature data sets tested so far, the temperature recovery law 
is: 

T = 8t x 

where 

T - extremum temperature minimum after a precipitation event 
t - recovery time in days 
x - exponent ranging between 0.62 < x < 0.92 


SIGNIFICANCE 


The correlations between diurnal temperature extrema for air and soil in 
drought and non-drought periods appear to follow different characteristic 
patterns, allowing an inference of soil moisture content from temperature 
data. The recovery of temperature extrema after a precipitation event also 
follows a characteristic exponential rise between two' limiting values which 
is an indicator of evaporation rates. If these indicators can be applied 
universally to regional temperature data, soil moisture content or drought 
conditions can be inferred directly from temperature measurements. 

FUTURE EMPHASIS 

An adequate statistical correlation testing of the 26 annual station data 
sets which have been assembled and plotted will be completed. All 26 data 

sets have already shown a similar characteristic temperature recovery after 
a precipitation event. HCMM data will be used for regional temperature data 
comparisons. An attempt will be made to incorporate this technique in esti- 
mates of ground water excesses related to seasonal variations of earth's 

rotation/L.O.D. estimates. This work will continue to be developed in con- 

sultation with Professor Helmut Landsfterg of tne University of Maryland. 
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Figure 1. Soil, Pan, ar.d Air Temperature with Precipitation and Evaporation for Titton, Georgia, 1979 













OF PCX 


nvh via id * ni dVA3 ivNama 


( «'pu|) NOIlVl'J03dd 


(d 0 ) 3an±vb3dW3± nos 
ivNama wrwixw; 


(do) NVd NOIlVdOdVA3 NI 3dniVb 3dl/\l31 
H31VM ONV biv IVNdma WniAIIXVW 


MKW * 

gL' 



t * * 4 c 




' • jfl - 


Figure 2. Soil, Pan, and Air Temperature with Precipitation and Evaporation for Tifton, Georgia. 1980 
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